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ABSTRACT 

Present day cosmic microwave background (CMB) studies require more accurate removal of Galactic 
foreground emission. This removal becomes even more essential for CMB polarization measurements. 
In this paper, we consider a way of filtering out the diffuse Galactic fluctuations on the basis of 
their statistical properties, namely, the power-law spectra of fluctuations. We focus on the statistical 
properties of two major Galactic foregrounds that arise from magnetized turbulence, namely, diffuse 
synchrotron emission and thermal emission from dust and describe how their power laws change 
with the Galactic latitude. We attribute this change to the change of the geometry of the emission 
region and claim that the universality of the turbulence spectrum provides a new way of removing 
Galactic foregrounds. For the Galactic synchrotron emission, we mainly focus on the geometry of the 
synchrotron emitting regions, which will provide useful information for future polarized synchrotron 
emission studies. Our model calculation suggests that either a one-component extended halo model or 
a two-component model, an extended halo component (scale height > Ikpc) plus a local component, 
can explain the observed angular spectrum of the synchrotron emission. For thermal emission from 
Galactic dust, we discuss general properties of a publicly available 94GHz total dust emission map 
and explain how we can obtain a polarized dust emission map. Based on a simple model calculation, 
we obtain the angular spectrum of the polarized dust emission. Our model calculation suggests that 
Ci oc Z^^^/^ for / > 1000 and a shallower spectrum for I < 1000. We discuss and demonstrate 
how we can make use of our findings to remove Galactic foregrounds using a template of spatial 
fluctuations. In particular, we consider examples of spatial filtering of a foreground at small scales, 
when the separation into CMB signal and foregrounds is done at larger scales. We demonstrate that 
the new technique of spatial filtering of foregrounds may be promising for recovering the CMB signal 
in a situation when foregrounds are known at a scale different from the one under study. It can also 
improve filtering by combining measurements obtained at different scales. 

Subject headings: MHD — turbulence — ISM:general — cosmic microwave background — Galaxy: struc- 
ture 



1. INTRODUCTION 

An important problem in the studies of the early uni- 
verse with the CMB fluctuations is related to separating 
them from Galactic foregrounds. The techniques of re- 
moving foregrounds are rather elaborate, but in most 
cases they include using the frequency templates of fore- 
ground emission. This requires multi-frequency measure- 
ments, which are not always available. Moreover, some 
foregrounds, e.g. the so-called spinning dust (Draine & 
Lazarian 1998ab; Finkbeiner et al. 2002; Lazarian & 
Finkbeiner 2003), demonstrate rather complex frequency 
dependence. Due to the utmost importance of obtaining 
CMB signal free of contamination, it is essential to con- 
sider other ways to remove foregrounds. One way to do 
this is to take into account the known spatial properties 
of emission. 

If a foreground has well-defined statistics of spatial 
fluctuations, one can devise techniques of removing the 
contribution of the foreground to the measured mi- 
crowave signal (see §2). The issue in this case is whether 
foregrounds have well defined behavior in terms of their 
spatial statistics. Determining this with the available 
data is the first thrust of our present study which we 
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pursue using the Galactic synchrotron and the Galac- 
tic dust emission data. The second thrust is devising 
possible ways of removing the foregrounds using the self- 
similarity of Galactic turbulence that gives rise to the 
foreground fluctuations. 

This work continues our brief study in Cho & Lazar- 
ian (2002a, henceforth CL02) where we argued that the 
properties of Galactic foreground radiation can be ex- 
plained on accepting that the interstellar medium (ISM) 
that provides the fluctuations is turbulent and therefore 
the spatial fluctuations of foregrounds inherit power-law 
spectra of the underlying magnetohydrodynamic (MHD) 
turbulence. Since that work, better understanding of the 
properties of the MHD turbulence has been achieved. For 
instance, it became clear that the spectrum of density in 
compressible MHD turbulence can be substantially shal- 
lower than the Kolmogorov spectrum that we assumed 
in CL02 (see Beresnyak, Lazarian & Cho 2005). More- 
over, the search for alternative procedures of removing 
foregrounds became more essential with the attempts 
to measure the polarization of the CMB radiation, es- 
pecially the enigmatic i3-modes. The latter motivates 
choice of foregrounds that we deal with in the present 
paper. Synchrotron and dust emissions are the sources 
of the polarized contamination for the CMB polarization 
studies. 

Diffuse Galactic synchrotron emission is an important 
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polarized foreground source, the understanding of which 
is essential for CMB studies, especially, in the range of 
10-100 GHz. It is known that the observed spectra of syn- 
chrotron emission and synchrotron polarization (see de 
Oliveira-Costa & Tegmark 1999 and references therein) 
reveal a range of power laws. The polarization of syn- 
chrotron emission traces magnetic fields and is perpen- 
dicular to the plane-of-sky magnetic field direction. Since 
the Galactic synchrotron emissivity is roughly propor- 
tional to the magnetic energy density, angular spectrum 
of synchrotron emission refiects st atist ics of magnetic 
field fluctuations in the Galaxy (see §A.ll for discussions). 

Thermal emission from dust is also an important 
source of polarized foreground emission in the range of 
frequencies larger than 60 GHz. The emission gets po- 
larized due to grain alignment (see Lazarian 2007 for a 
review). Therefore, polarization of dust, similar to syn- 
chrotron polarization, traces magnetic held fluctuations. 

What is the cause of the magnetic held fluctuations? 
As magnetic field lines are twisted and bend by turbulent 
motions in the Galaxy it is natural to think of the tur- 
bulence as the origin of the magnetic field fiuctuations. 
In fact, several earlier studies addressed the relation be- 
tween turbulence and the diffuse synchrotron foreground 
radiation. Tegmark et al. (2000) suggested that the spec- 
tra may be relevant to Kolmogorov turbulence. Chep- 
urnov (1999) and CL02 used different approaches, but 
both showed that the angular spectrum of synchrotron 
emission reveals Kolmogorov spectrum (C; ex l~^^'^) for 
larg43 values of multipole I. However, they noted that 
the spectrum can be shallower than the Kolmogorov one 
for smaller values of multipole /, due to density stratifi- 
cation in the halo (Chepurnov 1999) or the Galactic disk 
geometry (CL02J3. 

Recent research of compressible magnetized turbulence 
suggests that the fluctuations may not be necessarily 
Kolmogorov, to start with (see Beresnyak et al. 2005; 
Kowal & Lazarian 2007). Nevertheless, both observa- 
tions and numerical studies conflrm the power-law de- 
pendence of turbulence, even if the spectrum differs from 
the Kolmogorov one. This is also supported by theory, 
which states the self-similarity of turbulence. It is this 
self-similarity that gives us hope for a successful removal 
of foregrounds arising from turbulent interstellar media. 

At the same time, it is known that the spatial spectral 
index of foregrounds measured for different Galactic lat- 
itudes may differ. If the reason for this is unknown, this 
may make the weeding out of foregrounds using spatial 
templates unreliable. CL02 identified these changes with 
the variations of the geometry of the emission region. 
Thus, for the same Galactic latitudes one should expect 
the same slope of the spatial spectrum of the fiuctuations, 
which, for instance, allows one to extend the power-law 
spectrum of fluctuations measured for low spherical har- 
monics to higher spherical harmonics. Potentially, if the 
geometry of the emission region is known, this allows us 

^ In homogeneous turbulence, 'large values of multipole' means 
I larger than Icr = t/0, where 8 (in radian) is the angular size of 
the farthest eddies. CL02, for example, discussed that Icr ~ 30 
for the Galactic halo. However, in inhomogeneous turbulence, the 
practical value for Icr is an order of magnitude larger than ~ 30 
(see Chepurnov 1999; CL02). 

* In both approaches larger emissivity towards the disk plane is 
employed. 



to predict the expected changes of the indesQ- In this pa- 
per we provide more support for the conjecture in GL02. 
In this paper, we first present general properties and 
structure functions of a publicly available synchrotron 
foreground emission map. Then we investigate what 
kinds of Galactic halo structures can produce the ob- 
served structure function (and therefore angular spec- 
trum), which will be useful for the study of polarized 
synchrotron foreground. We also present the properties 
of a publicly available model dust emission map. Dust 
emission is one of the most important sources of polarized 
foreground radiation. Therefore, measurement of angu- 
lar power spectrum of such foreground is of great inter- 
est. Thus, we provide estimation of angular spectrum of 
polarized emission by foreground dust. This result is of 
great importance in view of recent interest to the CMB 
polarization. In §2 we explain a way of spatial removal 
of foreground emission and provide the summary of the 
expected scaling of foreground fluctuations arising from 
Galactic turbulence. In §3, we present statistical anal- 
ysis of the Haslam map, which is dominated by diffuse 
Galactic synchrotron emission. In §4, we investigate po- 
larized emission from thermal dust. In §5, we discuss how 
to utilize our knowledge to remove Galactic foregrounds. 
We provide the discussion of our results in §6 and the 
summary in §7. In Appendix, we review a simple model 
of the angular spectrum of synchrotron emission arising 
from MHD turbulence. In Appendix, we also present 
calculations of high-order structure functions of the syn- 
chrotron and the dust maps and we compare the results 
with those of turbulence. 

2. MOTIVATION: A NEW TECHNIQUE OF FOREGROUND 
REMOVAL 

2.1. Spatial Removal of Foregrounds 

Let us illustrate a possible procedure of the removal of 
Galactic foregrounds from the CMB signal. The cosmic 
microwave signals consists of the CMB signal I'^^^b and 
foregrounds I^ . When we correlate the microwave signal 
at points "1" and "2" , we get 

{{I?''^+I(){l§''^+Ii)) = {I^^'^I^'"^) + {I[li)^ 

(1) 
where we assume /C^^^ and I^ are uncorrelated. There- 
fore, the measured angular spectrum (jjneasured jg j^g^- 

the sum of Cf^^ and Cf and we have 
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CMB 



r^7neasured 



cr. 



(2) 



That is, if we know the angular spectrum of foregrounds 
cf , we can obtain the CMB angular spectrum Cf^^^ . 

How can one obtain Cf ? The well tested way of doing 
this is to use multi-frequency measurements of the CMB 
-I- foreground emission and separate the two components 
using the frequency templates of foregrounds. This ap- 
proach requires many measurements at different frequen- 
cies. In addition, for some foregrounds the frequency 
templates may be difficult to obtain. The so-called " spin- 
ning dust" foreground introduced in Draine & Lazar- 
ian (1998ab) presents an example of such a difficult-to- 
remove foreground. We also mention that the measure- 
ments of foregrounds at different frequencies may have 

^ Inverting arguments in CL02 one can use the changes of the 
foreground spectra to model the geometry of the emitting volume. 
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different spatial resolutions and the use of the maps with 
different resolutions may also present a problem. 

In this paper we address a somewhat different problem, 
which in its extremqj can be formulated in the following 
way. Imagine that we separated the foreground and the 
CMB signals at low resolution liow using the traditional 
multifrequency approach. Is it possible to use this in- 
formation to remove the foreground contribution from 
(jmeasured jTq^. ^ ^ //q^? For instance, the measurements 
of the Wilkinson Microwave Anisotropy Probe (WMAP) 
provide a high accuracy measure of Cf over a limited 
range of scales. If we know Cj^ as a function of Uow at 
scales smaller than those measured by the WMAP, then 
one can extrapolate C[ to / > Uow ■ These values of Cf 
can be used to filter the microwave measurements from 
balloon-born experiments using the procedure given by 
Eq. (12). Note that the balloon-born experiments usually 
have higher spatial resolutions, but not enough frequency 
coverage to remove foregrounds using the frequency tem- 
plates. 

The key question is to what extend we can predict Cf 
over a range of scales which is different from the range 
of scales at which Cf was measured. The answer to this 
question is trivial if Cf is a simple power law. While 
the actual spectra of foregrounds are more complex, in 
this paper we provide both theoretical arguments and 
the analysis of the foreground data which support the 
notion that C[ can be successfully extended beyond the 
range of I which is measured. 

We should stress that the filtering above is different 
from the accepted techniques of foreground removal us- 
ing frequency-dependent templates. The outcome of the 
latter procedures are maps of the foreground radiation 
and the CMB radiation. The filtering described above 
is of statistical nature. The result of it is Cp*^^ rather 
than emission map^. 

In view of above, it is important to determine to what 
extend the spatial properties of C[ are predictable, in 
particular, explore to what extend the power-law approx- 
imation is applicable. This paper provides a study of spa- 
tial statistical properties of fluctuations of synchrotron 
and dust emission and relates those to the properties of 
the underlying turbulence. It also provides an example 
of the filtration procedure that we advocate. 

Another important issue is to determine the reasons 
for the change of the power law behavior with latitude. 
This is what we study below for the synchrotron and 
dust foreground emissions. The ultimate goal of this re- 
search is to obtain models of Galactic foregrounds which 
provide a good fit for the foreground spatial spectrum at 
arbitrary scales and arbitrary latitudes. This paper is a 
step to constructing such a model. 

2.2. Power-law behavior of interstellar turbulence 

For the spatial removal procedure to work, we should 
understand the spatial spectra of foregrounds. Since the 

® Less extreme cases would involve the use of the known spa- 
tial properties of Cf to increase the accuracy of the removal of 
foregrounds within traditional techniques. We do not discuss these 
more sophisticated procedures in this paper. 

"^ It is easy to see that the phase information of foreground emis- 
sion is lost in the process of such a filtering. However this infor- 
mation is not necessary for the Cp**^^ recovery. 



spatial fluctuations of foregrounds inherit the spectra of 
the underlying interstellar turbulence, we summarize the 
spectral behavior of interstellar turbulence. 

It is generally accepted that the interstellar medium 
(ISM) is magnetized and turbulent (see reviews by 
Elmegreen & Scalo 2004 and McKee & Ostriker 2007). 
The so-called "Big Power Law in the Sky" correspond- 
ing to the Kolmogorov-type turbulence was reported in 
Armstrong, Rickett, & Spangler (1995). Recently this 
law based on the measurements of the radio scintilla- 
tions arising from electron density inhomogeneities has 
been extended to larger scales using the WHAM H^ fluc- 
tuations (Chepurnov & Lazarian 2010). 

Spectra of magnetic turbulence has been studied us- 
ing Faraday rotation measurements (see Haverkorn et al. 

2008) and starlight polarization (Hildebrand et al. 2009). 
The interpretation of the results are more challenging in 
those cases. 

Molecular and atomic spectral lines present a very 
promising way of studying turbulence. The interstel- 
lar lines are known to be Doppler-broadened due to 
turbulent motions. Obtaining spectra from Doppler- 
broadened lines is not trivial, however. A lot of research 
in this direction based on the use of the so-called "ve- 
locity centroids" , which were the main tool to study 
velocity fluctuations, has been shown to produce erro- 
neous results for supersonic turbulence (Lazarian & Es- 
quivel 2003; Esquivel & Lazarian 2005; Esquivel et al. 
2007). At the same time new techniques based on the 
theoretical description of the Position-Position- Velocity 
(PPV) data cubes, namely, the Velocity Channel Analy- 
sis (VCA) and the Velocity Coordinate Spectrum (VCS) 
have been developed (Lazarian & Pogosyan 2000, 2004, 
2006, 2008), tested (see Stanimirovic & Lazarian 2001; 
Lazarian, Pogosyan, & Esquivel 2002; Esquivel et al. 
2003; Chepurnov & Lazarian 2009; Padoan et al. 2006, 

2009) and applied to the observational data to obtain 
the characteristics of the velocity turbulence (see Lazar- 
ian 2009 for a review). 

VCA and the VCS techniques reveal that the veloci- 
ties for interstellar turbulence may be somewhat steeper 
than the Kolmogorov one , while the density spectra of 
the fluctuations may be substantially more shallow than 
the spectrum of Kolmogorov fluctuations (see Padoan et 
al. 2006, 2009; Chepurnov et al. 2010). This agrees 
well with the numerical studies of the magnetized su- 
personic turbulence (Beresnyak, Lazarian & Cho 2005; 
Kowal, Lazarian & Beresnyak 2007; Kowal & Lazarian 
2010). For the subsonic turbulence the spectrum of mag- 
netized media gets the values close to the Kolmogorov in- 
dex (see Goldreich & Sridhar 1995; Lazarian & Vishniac 
1999; Cho & Vishniac 2000; MuUer & Biskamp 2000; 
Maron & Goldreich 2001; Lithwick & Goldreich 2001; 
Cho, Lazarian & Vishniac 2002; Cho & Lazarian 2002b, 
2003; Boldyrev 2006; Beresnyak & Lazarian 2006, 2009). 
In view of that we believe that the Kolmogorov spectrum 
can be used as a proxy of the underlying spectra of veloc- 
ity and magnetic fleld, while a more cautious approach 
should be demonstrated to density in highly compressible 
environments, e.g. molecular clouds. 

Below we shall show that the spectra of turbulence de- 
rived from the analysis of foreground is consistent with 
both the results of dedicated observations of turbulence 
as well as the theoretical expectation for the MHD tur- 
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2.3. Data Sets 



We use the 408MHz Haslam all-sky map (Haslam et 
al. 1982) and a model 94GHz dust emission map that 
are available on the NASA's LAMBDA website 0. Both 
maps were reprocessed for HEALPix (Gorski et al. 2005) 
with nside=512 (7' resolution). 

The original Haslam data were produced by merging 
several different data-sets. "The original data were pro- 
cessed in both the Fourier and spatial domains to miti- 
gate baseline striping and strong point sources" (see the 
website for details) . The angular resolution of the orig- 
inal Haslam map is ^ 1°. Galactic diffuse synchrotron 
emission is the dominant source of emission at 408MHz. 

The 94 GHz dust emission map is based on the work 
of Schlegel, Finkbeiner, & Davis (1998) and Finkbeiner, 
Davis, & Schlegel (1999). Schlegel et al. (1998) com- 
bined 100/i?7i maps of IRAS (Infrared Astronomy Satel- 
lite) and DIRBE (Diffuse Infrared Background Exper- 
iment on board the COBE satellite) and removed the 
zodiacal foreground and point sources to construct a full- 
sky map. Finkbeiner et al. (1999) extrapolated the 
lOOfim emission map and 100/240/i77i flux ratio maps 
to sub-millimeter and microwave wavelengths. The 94 
GHz dust map we used is identical to the two-component 
model 8 of Finkbeiner et al. (1999). The angular reso- 
lution of the 100-micron map is ~ 6' and that of the 
temperature correction derived from the 100/240-micron 
ratio map is r^ 1° (Finkbeiner et al. 1999), which corre- 
sponds to I - 180761° - 180. 

3. SPATIAL STATISTICS OF DIFFUSE GALACTIC 
SYNCHROTRON EMISSION 

In this section, we analyze the Haslam 408MHz all-sky 
map, which is dominated by Galactic diffuse synchrotron 
emission. Our main goal is to explain the observed syn- 
chrotron angular spectrum using simple turbulence mod- 
els. The result in this section will be useful for sophisti- 
cated modeling of polarized synchrotron emission. 

3.1. General properties of diffuse Galactic synchrotron 

emission 

In this section, we study synchrotron emission from the 
Galactic halo (i.e. b > 30°) and the Galactic disk (i.e. 
I&l < 2°) separately. Our main goal is to see if statistics of 
synchrotron emission from the halo is consistent with tur- 
bulence models. When it comes to synchrotron emission 
from the Galactic disk, it is not easy to separate diffuse 
emission and emission from discrete sources. Therefore, 
we do not try to study turbulence in the Galactic disk. 
Instead, we will try to estimate which kind of emission 
is dominant in the Galactic disk. 

There exist several models for the diffuse Galactic ra- 
dio emission. Beuermann, Kanbach, & Berkhuijen (1985) 
showed that a two-component model, a thin disk embed- 
ded in a thick disk, can explain observed synchrotron 
latitude profile. They claimed that the equivalent width 
of the thick disk is about several kiloparsecs and thin disk 
has approximately the same equivalent width as the gas 
disk. They assumed that, in the direction perpendicular 

® |http://lambda.gsfc.nasa.gov/| 



to the Galactic plane, the emissivity e of each component 
follows 

e(z) = e{0)sech\z/zo), (3) 

where z is the distance from the Galactic plane and e(0), 
b, and zq are constants. The half-equivalent width of the 
disk, which is proportional to zq, at the location of the 
Sun is ~ 2kpc. 

Recently, several Galactic synchrotron emission mod- 
els have been proposed in an effort to separate Galac- 
tic components from the WMAP polarization data (see, 
for example. Page et al. 2007; Sun et al. 2008; Miville- 
Deschenes et al. 2008; Waelkens et al. 2008). All the 
models mentioned above assume the existence of a thick 
disk component with a scale height equal to 1 kpc. Sun et 
al. (2008) considered an additional local spherical com- 
ponent motivated by the local excess of the synchrotron 
emission that might be related to the "local bubble" (see, 
for example, Fuchs et al. 2008). 

The detailed modeling of the Galactic synchrotron 
emission is beyond the scope of our paper. We will sim- 
ply assume that there is a thick component with a scale 
height of ^ Ikpc. We will also assume that there could 
be an additional local spherical component. Then, in 
the following subsections, we will consider the relation 
between spectrum of 3-dimensional turbulence and the 
observed angular spectrum of synchrotron emission. 

3.2. Structure function of the ^08MHz Haslam map 

In Appendix, we discussed the relation between the 3D 
spatial MHD turbulence spectrum and the observ ed 2D 
angular spectrum of synchrotron emission (see Eq. |All| 
for a quick summary). But, in Appendix we assumed 
that the emission arises from a spherical region filled with 
homogeneous turbulence. In this section, we will show 
that the modulation of synchrotron intensity of the emit- 
ting volume can also affect the observed angular spec- 
trum of synchrotron emission. 

Earlier studies showed that the angular spectrum 408- 
MHz Haslam map has a slope close to —3: Ci oc l~^ 
(Tegmark & Efstathiou 1996; Bouchet, Gispert, & Puget 
1996). Recently La Porta et al. (2008) performed a com- 
prehensive angular power spectrum analysis of all-sky 
total intensity maps at 408MHz and 1420MHz. They 
found that the slope is close to -3 for high Galactic lati- 
tude regions. Other results also show slopes close to -3. 
For example, using Rhodes/HartRAO data at 2326 MHz 
(Jonas, Baart, & Nicolson 1998), Giardino et al. (2001b) 
obtained a slope ^ 2.92 for high Galactic latitude re- 
gions with |6| > 20°. Giardino et al. (2001a) obtained 
a slope ~ 3.15 for high Galactic latitude regions with 
l&l > 20° from the Reich & Reich (1986) survey at 1420 
MHz. Bouchet & Gispert (1999) also obtained a slope 
~ l~^ spectrum from the 1420 MHz map. 

In general, synchrotron emission from the Galactic disk 
makes it difficult to measure the angular spectrum of 
synchrotron emission from the Galactic halo. In order to 
avoid the contamination by the Galactic disk, one may 
mask out the low Galactic latitude regions. This can be 
done, for example, by setting all synchrotron intensity 
to zero for regions with \b\ < hcut- However, the angu- 
lar spectrum obtained with the Galactic mask exhibits 
spurious oscillations. Moreover, the spectrum obtained 
with a mask may not be the true one because it is con- 
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Fig. 1. — Haslam 408MHz map. Left: Second-order angular structure function (for b > 30°) shows a slope of ~ 1.2, which is shallower 
than that of Kolmogorov turbulence (slope = 5/3). Right. Structure function as a function of Galactic latitude. From bottom to top, the 
second-order angular structure functions are obtained for thin stripes (|A6| < 2°) along Galactic latitudes of 30°, 10°, and 0°. 



taminated by the mask. In principle, one may correct 
such oscillations and estimate the true spectrum, using 
the convolution theorem: the Fourier coefficients (or, in 
this case, spherical harmonic coefficients) of the masked 
data are convolution of those of the true data and those 
of the mask. However, the practical implementation of 
the method is not simple. 

Another approach is to estimate 2-point correlation 
function first and to extract the angular power spectrum 
from it (Szapudi et al. 2001; see also Eq. [in])- This 
method is free of artifacts caused by the mask. But the 
angular spectrum Ci obtained in this way is, in general, 
noisy and it requires a lot of calculations to accurately 
measure the slope of the spectrum. 

We are interested in the slope of the angular spectrum 
on small angular scales, which is the same as that of 
the underlying 3D spatial turbulence spectrum (see Ap- 
pendix). Therefore, in this paper, we use yet another 
approach. We first calculate the second-order angular 
structure function: 



i?2(e)=<|/(ei)-/(e2)|2>. 



(4) 



where /(e) is the intensity of synchrotron emission, ei 
and 62 are unit vectors along the lines of sight, 9 is the 
angle between ei and e2, and the angle brackets denote 
average taken over the observed region. Then, we extract 
the slope of the angular spectrum using the relation 

£'2(6') oc 6l™-2 ^ Q oc r™ ^ £^315 oc fc-" (5) 

for small angular scales (see ijA.2l) . 

We note that excessive care is required in the presence 
of white noise. In the presence of white noise, second- 
order structure function will be D2{9) =< |/(ei) -\- 5i — 
/(e2) - <52|' >-< |/(ei) - /(e2)|2 > + < \6i ~ 52\'' > 
, where 5i and 82 represent noise (Chepurnov, private 
communication). If the second term on the right (< |(5i — 
(52 P >) is sufficiently smaller than the first term on the 
right (< |/(ei) — /(e2)|^ >), we can ignore the noise. 
Otherwise, the noise can interfere accurate measurement 
of the slope. In the Haslam map, the level of the white 
noise seems to be negligible. The reason is as follows. 
Our measurements show that 1)2(0.015°) ~ 0.05 in the 
Haslam map. This means that the < |5i — ^2^ > term 
is no larger than 0.05, which is sufficiently smaller than 



values D2{9) shown in Fig. [T] 

In the left panel of Fig. [1] we show the second-order 
structure function for the Galactic halo (i.e. |6| > 30°). 
The slope of the second-order structure function lies be- 
tween those of two straight lines. The steeper line has a 
slope of 4/3 and the other one has a slope of 1. The actual 
measured slope is ^ 1.2. This result implies that the 3D 
turbulence spectrum is E^jy^k) oc k~^'^, which is shal- 
lower than the Kolmogorov spectrum E^oik) oc k~^^^^. 

Now a question arises: why is the slope shallower than 
that of Kolmogorov? Chepurnov (1999) provided a dis- 
cussion of the effects of density stratification on the slope. 
He used a Gaussian disk model and semi-analytically 
showed that the slope of the angular spectrum can be 
shallower than that of Kolmogorov. In the next subsec- 
tion, we present our understanding of the effect. 
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Fig. 2. — Halo models with stratification. Upper plot: Ex- 
ponentially stratified halo. We take zo(= rg) = Ikpc and L 
(=eddy size) = lOOpc. Lower plot: Two-component halo. We take 
21 (= ri) = Ikpc, Z2{= ^2) = lOOpc, and L (=eddy size) = lOOpc. 



3.3. Model calculations 

Below we take into account that the emission does 
depend on the distance from the Galactic plane. Syn- 
chrotron emission models (see previous subsection) as- 
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sume either exponential (e ^/^o) or square of hyperbolic 

secant (sech [z/zq]) l^w for synchrotron emissivity, where 
z is the distance from the Galactic plane and 6 is a con- 
stant. For simplicity, we assume the observer is at the 
center of a spherical halo. That is, the geometry is not 
plane-parallel, but spherical. In what follows, we use r, 
instead of z, to denote the distance to a point. 

To illustrate the effects of this inhomogeneity, we test 
3 models: 

1. Homogeneous halo: Turbulence in halo, thus 
emissivity, is homogeneous. Turbulence has a sharp 
boundary at dmax = Ifcpc. The outer scale of tur- 
bulence is lOOpc. Basically, this model is the same 
as the one we considered in ijA) 

2. Exponentially stratified halo: Emissivity 
shows an exponential decrease, e(r) oc e"'"''"''. We 
assume rg = Ifcpc and the halo truncates at r = 
8fcpc. The outer scale of turbulence is lOOpc. See 
Fig.H 

3. Two-component halo: Emissivity decrease as 



where €2 



lOei, 



e(r) oc eie"''/''! -I- e2e"''/'■^ 
ri = Ifcpc, and r2 = lOOpc. The halo truncates at 
r = 8fcpc. The outer scale of turbulence is lOOpc. 
The second component mimics local enhancement 
of synchrotron emissivity. See Fig. [2] 

We numerically calculate the angular correlation func- 
tion w(G) and the second-order structure D2{9) from 



w{e) = I dn I dr2 /C(|ri - r2|)e(ri)e(r2), (6) 
D2{e)^T-w{e), (7) 



where |ri — r2| = r^ 



2rir2COS0, e(r) is the syn- 



chrotron emissivity, T = Y\uYg^Qw{9), and wc use the 
spatial correlation function 1C{r) obtained from the rela- 
tion: 

IC{r) <x r 47Tk^E3D{k)^^^dk (8) 

Jo kr 

where the spatial spectrum of emissivity Eio has the 
form: 

, , f constant ii k < ko , , 

which is the same as Kolmogorov spectrum for k > ko {'^ 
1/L). The reason we use a constant spectrum for k < ko 
is explained in Appendix B (see also Chepurnov 1999). 
We obtain the angular spectrum from the relation: 



Ci oc Pi{cose)K{cose) d{cosd), 



(10) 



where Pi is the Legendre polynomial. 

In Fig. [21 we plot the calculation results. The angular 
correlation function w{0) does not change much when 
9 is small, and follows ~ (tt — 9)/sm6 ^ 1/6 when 9 
is large. The critical angle is a few degrees for homo- 
geneous model (thick solid curve) and single-component 
exponential model (dotted curve). As we discussed ear- 
lier, the critical angle for homogeneous turbulence is 
^ {L/dmaxY'"'' ~ 6°, where d^ax (= ^kpc in our model) 
is the distance to the farthest eddy. In Fig. |3] (left panel) 



we clearly see that the slope of w{9) changes near 6* ~ 6°. 
The second-order structure function D2{9) also shows a 
change of slope near the same critical angle {9 ~ 6°). 
In single-component exponential model (dotted curve), 
the value of dmax is not important. Instead, the scale 
height zq is a more important quantity, which is \kpc in 
our model. In left and middle panels of Fig. [31 we ob- 
serve that the single-component exponential model also 
show a change of slope near 9 ^ a. few degrees. There- 
fore, we can interpret that the critical angle for stratified 
turbulence is '^ L/ zq, instead of ^^ L/d^ax 

Now, it is time to answer our earlier question of why 
the observed slope is shallower than that of Kolmogorov. 
Let us take a look at the right panel of Fig. [3] All 3 
models show that the slope of Ci is almost Kolmogorov 
one for I > 200 > Icr '^ T^dmax/ L ^ 30. However, if we 
measure average slope of q between I — 10 and 200, we 
obtain slopes shallower than Kolmogorov. The single- 
component model and the homogeneous model give sim- 
ilar average slopes around —3. However, the homoge- 
neous model gives a more abrupt change of slope near 
/ ^ 30. In fact, right panel of Fig. [31 shows a notice- 
able break near ^ ^ 50. The average slope of the two- 
component model gives more or less gradual change of 
the slope and the observed slope is very close to —3 for 
a broad range of multipoles I. It is difficult to tell which 
model is better because the models are highly simpli- 
fied. Nevertheless, the two-component model looks the 
most promising, which is not so surprising because two- 
component model has more degree-of-freedom. 

Note that, compared with the two-component model, 
the single-component model shows a more or less sudden 
change of slope near I ^ Icr ^ t^zq/ L ~ 30. Therefore, if 
the single-component model is correct, the scale height zq 
cannot be much larger than '^ 10 times the outer scale 
of turbulence L. If zq is much larger than ~ lOL, Icr 
becomes smaller and we will have almost Kolmogorov 
slope for I > 10. We also note that it is possible that 
3D spatial turbulence spectrum itself can be shallower 
than the Kolmogorov one. That is, it is possible that 
spectrum of B(r), hence that of B^{r), can be shallower 
than the Kolmogorov one. For example, some recent 
studies show that strong MHD turbulence can have a ^ 
^-3.5 spectrum, rather than k~^^'^ (Maron & Goldrcich 
2001, Boldyrev 2005; Beresnyak & Lazarian 2006JB- If 
this is the case, the observed angular spectrum can be 
slightly shallower than the Kolmogorov for I > Icr- 

To summarize this subsection, our simple model calcu- 
lations imply that C/ will be very close to the underlying 
spectrum of magnetic turbulence for large values oil {I > 
a few time 100). The corresponding spectral slope is ex- 
pected to be close to the Kolmogorov one (see §3). For 
intermediate values of / (e.g. 10 < ^ < 200), the aver- 
age slope is shallower than the Kolmogorov one. Thus 
our modeling shows the consistency of the observational 
spectra with the expectations. Studies of the observed 
spectra at higher I may be useful for better testing of our 
predictions. 



^ The reason for the spectrum being shallow in simulations is 
unclear. It may also be the result of the limited dynamical range 
in the presence of non-locality of MHD turbulence (Beresnyak & 
Lazarian 2009). 
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Fig. 3. — Model calculations. Three toy models for emissivity profiles in the Galactic halo are considered: homogeneous (solid line) 
halo, exponentially stratified halo (dotted line), and two-component exponential halo (dashed line). Left: Angular correlation functions. 
When angular separation is large, angular correlation functions follow the universal relation: (tt — 9)/sin9 ~ 1/6. Middle: Second-order 
structure functions. When angular separation is small, the slope of the homogeneous turbulence (solid line) is compatible with 5/3. But 
those of stratified halo models are shallower. The two-component model (dashed line) shows a slope compatible with —3, if we measure 
average slope between 8 ~ 0.5° and ~ 10°. Right: Angular spectra. All spectra are compatible with the 3D spatial turbulence spectrum 
of l~-'-^'^ for I > a few time 100. The stratified halo models show shallower slopes, if we naively attempt to fit it with a single power law 
and measure the resulting averaged slope between Z ~ 10 and ~ 200. The homogeneous turbulence model also gives a shallower slope for 
these values of I. But, its spectrum shows a distinct break near I ~ 50. 



3.4. Synchrotron emission from Galactic disk 

In right panel of Fig. [T] we show how the second-order 
structure function changes with Galactic latitude. The 
lower curve is the second-order angular structure func- 
tion obtained from pixels in the range of 28° < b < 32°. 
The middle and upper curves are the second-order angu- 
lar structure functions obtained from pixels in the range 
of 8° < 6 < 12° and -2° < 6 < 2°, respectively. The 
middle and upper curves clearly show break of slopes 
near ~ 3° and ~ 1.5°, respectively. When the angu- 
lar separation is larger than the angle of the break, the 
structure function becomes almost flat. As we move to- 
wards the Galactic plane, the sudden changes of slopes 
happen at smaller angles. 

There are at least two possible causes for the break 
of slope. First, a geometric effect can cause it. As we 
discussed in the Appendix, change of slope occurs near 
9c '^ L/dmax- As we move towards the Galactic plane, 
the distance to the farthest eddy, dmax, will increase. As 
a result, the critical angle 9c ^ L/dmax will decrease. 
Therefore, we will have smaller 9c towards the Galac- 
tic plane. This may be what we observe in the right 
panel of Fig. [TJ Second, discrete synchrotron sources 
can cause flattening of the structure function on an- 
gular scales larger than their sizes. Although the map 
we use was reprocessed to remove strong point sources, 
there might be unremoved discrete sources. When fila- 
mentary discrete sources dominate synchrotron emission, 
the second-order structure function will be flat on scales 
larger than the typical width of the sources. In reality, 
both effects may work together. 

In view of the variations of the spatial spectral slope 
of the synchrotron emission at low Galactic latitudes, 
the use of the foreground removal procedure discussed in 
§2 is more challenging for those latitudes. At the same 
time, the high-Z fluctuations corresponding high latitudes 
should be possible to remove reliably with the procedure 
in §2 due to the observed regular power-law behavior. 

3.5. On the polarized synchrotron emission 

Roughly speaking, the shape of the angular spectrum 
of polarized synchrotron emission will be similar to that 
of the total intensity at mm wavelengths. However, it 
is expected that at longer wavelengths, Faraday rotation 



and depolarization effects should cause flattening of the 
angular spectrum, which has been actually reported (see 
de Oliveira-Costa et al. 2003 and references therein). 
On the other hand. La Porta et al. (2006) analyzed the 
new DRAO 1.4GHz polarization survey and obtained an- 
gular power spectra with power-law slopes in the range 
[—3.0,-2.5]. More observations on the polarized syn- 
chrotron foreground emission can be found in Ponthieu 
et al. (2005), Giardino et al. (2002), Tucci et al. (2002), 
Baccigalupi et al. (2001). In this paper, we do not dis- 
cuss the properties of the polarized synchrotron emission. 
Readers may refer to recent models about the polarized 
synchrotron emission (Page et al. 2007; Sun et al. 2008; 
Miville-Deschenes et al. 2008; Waelkens et al. 2008). 

4. POLARIZED EMISSION FROM DUST 

Polarized radiation from dust is an important com- 
ponent of Galactic foreground that strongly interferes 
with the intended CMB polarization measurements (see 
Lazarian & Prunet 2001). Therefore, the angular spec- 
trum of the polarized radiation from the foreground dust 
is of great interest. One of the possible ways to estimate 
the polarized dust radiation at the microwave range is to 
measure star-light polarization and use the standard for- 
mulae (see, for example, Hildebrand et al. 1999) relating 
polarization at different wavelengths. This approach in- 
volves a number of assumptions the accuracy of which we 
analyze below. In this section, we describe how we can 
obtain a map of polarized dust emission using starlight 
polarization and we discuss the angular spectrum of the 
polarized foreground emission from thermal dust at the 
high Galactic latitude (say, |6| > 20°). 

4.1. Properties of the 94GHz Dust Emission Map 

Let us begin with a model dust emission map 
created by Finkbeiner et al. (1999), which 

is available at the NASA LAMBDA website 
(|http://lambda.gsfc.nasa.gov/). As we will explain 
later in this section, we can derive a polarized intensity 
map from this kind of total intensity map. 

We note that the difference between the dust emission 
and synchrotron emission (discussed in the previous sec- 
tion) is expected. The origin of synchrotron emission 
is related to cosmic ray electrons which are distributed 
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within an extended magnetic halo (Ginzburg & Ptuskin 
1976). At the same time, dust is expected to be locahzed 
mostly within the Galactic plane. In Fig. |4l we present 
statistical properties of the map. The map shows rough 
constancy of emission for high Galactic latitude region 
when multiplied by sin b (left panel of Fig. HJ . The sin b 
factor also appears in the PDF: the sin 6 factor makes the 
PDF more symmetric (right panel of Fig.|4]). Therefore it 
is natural to conclude that a disk component dominates 
the dust map. 

As in the Haslani map, we do not try to obtain angular 
spectrum of dust emission map directly. Instead, we use 
the second-order structure function in order to reveal the 
angular spectrum on small angular scales. Indeed, the 
second-order structure function of the dust map (see [JdJ 
shows a slope of ~ 0.6, which corresponds to angular 
spectrum of ~ l~'^-^. 

The slope of the angular power spectrum of the model 
dust emission map is very similar to that of the original 
FIR data. Schlegel et al. (1998) found a slope of -2.5 
for the original FIR data. On the other hand, other 
researchers found slopes close to —3 from other observa- 
tions (see Tegmark et al. 2000 and references therein; see 
also Masi et al. 2001). 

Dust density fluctuations mostly arise from cold dense 
phases of interstellar medium (see Draine & Lazarian 
1998 for a list of idealized ISM phases). There the turbu- 
lence is known to be supersonic, which is vividly revealed, 
for instance, by Doppler broadening of observed molecu- 
lar lines from Giant Molecular Clouds (GMCs) (see Mc- 
Kee & Ostriker 2007). The shallow spectrum of den- 
sity fluctuations is consistent with the numerical simula- 
tions of supersonic MHD turbulence in Beresnyak et al. 
(2005). The shallow spectrum was later also reported in 
supersonic hydro turbulence in Kim & Ryu (2005), which 
indicates that the effect is not radically changed by the 
magnetic field. The latter makes the conclusion about 
the shallow spectrum independent on the degree of the 
interstellar medium magnetization and sub-Alfvenic ver- 
sus super-Alfvenic character of turbulence there. Thus 
we claims the observed spectra should be associated with 
the shallow spectra of underlying density fluctuations in 
the denser part of the interstellar medium. We predict 
that the shallow power law arising from density fluctu- 
ations extends from the scales of the turbulence energy 
injection to the dissipation scales. As a result, power-law 
extension of the observed data and the corresponding fil- 
tering using Eq. ([2]) is possible. 

4.2. Map of polarized dust emission from starlight 
polarization 

In general, it is advantageous to use all possible sources 
of information about foregrounds in order to improve 
their removal. In this section we discuss how the starlight 
polarization maps can be used to construct the maps 
of dust polarized emission. In principle, we can con- 
struct a polarized dust emission map at mm wavelengths 
{Ipoi,mmil, b)) from a dust total emission map {Immil, b)) 
and a degree-of-polarization map {Pem,mm(J,b)) at mm 
wavelengths: 

where {l,b) denotes the Galactic coordinate. However, 
neither Ipoi,mm{hb) nor Imm{lTb) is directly available. 



Therefore, we need indirect methods to get Ipoi.mm{l,b) 
and Immil, b). 

Obtaining a dust total emission map {Imm{l, b)) is rel- 
atively easy because dust total emission maps at FIR 
wavelengths are already available from the IRAS and 
COBE/DIRBE observations. Using the relation 



i{l,b) = Iwof^m{l,b){lmm/lQQ^im) ^, 



(12) 



where 1^/3^2, one can easily obtain an emission map 
at mm wavelengths (Imm) from the maps at 100 fim or 
240 fxm. However, more sophisticated model dust emis- 
sion maps at mm wavelengths already exist. For exam- 
ple, Finkbeiner et al. (1999) presented predicted full-sky 
maps of microwave emission from the diffuse interstellar 
dust using FIR emissions maps generated by Schlegel et 
al. (1998). In fact, the model dust emission map we an- 
alyzed in the previous section ( §4.1|) is one of the maps 
presented in Finkbeiner et al. (1999). Therefore, we can 
assume that the thermal dust emission map {Imm{l,b)) 
is already available. 

Obtaining a degree-of-polarization map at mm wave- 
lengths {Pem,mm(jTb)) is relatively more complicated. 
We can use measurements of starlight polarization at 
optical wavelengths to get Pem,7nm{l,b). The basic 
idea is that the degree of polarization by emission at 
mm {Pem,mm) IS related to that at optical wavelengths 
(Pem, optical), which in tum is related to the degree 
of polarization by absorption at optical wavelengths 

\-^ahs , optical ) ' 

^ahs^optical ^ ^cm^optical ^ ^em.mm- V *^j 

We describe the relations in detail below. 

When the optical depth is small, we have the following 
relation (see, for example, Hildebrand et al. 2000): 



P. 



ein,opt 



-t'^abSjOpt/ "^1 



(14) 



where P. 



em,opt 



is the degree of polarization by emis- 
sion and T is the optical depth (at optical wavelengths). 
We obtain polarization by emission at mm wavelengths 
{Pem.mm) using the relation 



^em,rnm — ^em,opt 



L-rj 



G„ 



L-r? 



C„ 



I 



Or, 



o„ 



o„ 



c„ 



(15) 
where C's are cross sections (of grains as projected on 
the sky) that depend on the geometrical shape (see, for 
example, the discussion in Hildebrand et al. 1999; see 
also Draine & Lee 1984) and dielectric function t — e\ + 
ie2 (see Draine 1985) of grains. 

For mm wavelengths, it is easy to calculate the ratio in 
Eq. (ITS]) because the wavelength A is much greater than 
the grain size a (i.e. A 3> 27ra). In this case, if grains are 
oblate spheroids with ai < 02 = 03 and short axes (ai) 
of grains are perfectly aligned in the plane of the sky, we 
have 



opt 



C, 



2ttV 



e2(A) 



where L values are defined by 

ii = [(l + /')//'][! -(l//)arctan/], 
i2=i3 = (l-ii)/2. 



(16) 



f = ia2/aiy 



1 



(17) 



Synchrotron and Dust Foregrounds 



Dust * sin(b' 



0.0100 



cn 

D 



> 0.0010 



0.000' 




20 



40 



60 



80 





2.5 




2.0 




1.5 


f-) 




CL 


1.0 




0.5 




0.0 



Dust (b>30°) 



ntensity 
ntensity*sin(b) 




0.001 0.010 

ntensity (mK) 



Fig. 4. — Average dust emission intensity times sin 6. The 94GHz dust emission map has been used. Left: Latitude profile of averaged 
dust emission intensity time sin b. Average value is roughly constant in high latitude region. The solid line represents the average taken over 
360°. The dashed line and the dotted line depict average taken near the Galactic center (—45° < loal < 45°) and the Galactic anti-center 
(135° < Ical ^ 225°), respectively. Right: PDF for |6| > 30°. PDF of dust intensity times sinb (dotted curve) shows a rough symmetry, 
but that of dust intensity (solid curve) does not. 
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Fig. 5. — The ratio of Pem, mm /Pem, optical- Left: The polarization ratio vs. the axis ratio (02/01) of aligned oblate spheroidal grains. 
The polarization ratio shows only a weak dependence on the axis ratio. We assume that the grain are perfectly aligned with their long 
axes perpendicular to magnetic field, the grain size is O.l/^m, magnetic field is perpendicular to the line-of-sight, XopUcal = 0.5/im, and 
Amm = lOOOfira. Right: The polarization ratio vs. the angle between magnetic field and the plane of the sky. We assume that the grain 
size is 0.1/im, the grain axis ratio is 1.5, ^optical = 0-5fj,m, and Xmm = lOOOfira. 



(see, for example, Hildebrand et al. 1999). 

However, for optical wavelengths, the condition A ^ 
27ra is not always valid and, therefore, the expression 
in Eq. ()16|) returns only approximate values. For accu- 
rate evaluation of the cross sections, one should use nu- 
merical methods. Fortunately, several numerical codes 
are publicly available for such calculations (for exam- 
ple, DDSCAT package by Draine & Flatau 1994, 2008; 
ampld.lp.f by Mishchenko 2000). We use ampld.lp.f 
to calculate the ratio in Eq. ([T5|) . We assume that 
the grains are oblate spheroids, grain size is 0.1/im, 
^optical — 0.5/im, and Xmm = 1000//m. Left panel of 
Fig. [S] shows that the ratio is around 1.5 when magnetic 
field is perpendicular to the line-of-sight. It also shows 
that the ratio of Pem.mm/ Pem, optical is almost indepen- 



dent of the grain axis ratio. 

In this subsection, we described a simple way to ob- 
tain a polarized map at mm wavelengths. However, ac- 
tual implementation of the method can be more compli- 
cated due to the following reasons. First, we used an 
assumption that the grains that produce optical absorp- 
tion produce also microwave emission. But, this is not 
true in general (see Whittet et al. 2008). Second, the ex- 
pressions in Eqs. (|15p and (J16p are valid when magnetic 
field direction is fixed and perpendicular to the line-of- 
sight and all grains are perfectly aligned with the mag- 
netic field. If this is not the case, Eq. (|15p will become 
Pem,mm oc Pem,opticai with the constant of proportion- 
ality that depends on magnetic field structure and the 
degree of grain alignment. The effect of partial align- 
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ment is expected to be less import anlF^. 

The effect of non-perpendicular magnetic field can be 
potentially important. We perform a numerical calcula- 
tion using ampld.lp.f to evaluate the effect. We assume 
that the grains are oblate spheroids, grain size is 0.1/im, 
^optical = 0.5/xm, and Xmm = 1000/^m. Right panel of 
Fig. [S] shows that the polarization ratio drops from ~ 1.5 
to ^ 1.1 when the angle (between magnetic field and the 
plane of the sky) changes from 90° to ^ 5°. Therefore, 
the effect is not very strong and can be potentially cor- 
rected foiEl. 

4.3. Angular spectrum of polarized emission from 
thermal dust 

After wc have constructed a map of the polarized emis- 
sion from thermal dust, we can obtain the angular spec- 
trum. However, if we are interested in only the shape 
of the angular spectrum, we do not need to construct 
the polarized thermal dust emission map. We can get 
the shape of the the angular spectrum directly from the 
starhght polarization map Pabs,opUcaiil,b). 



Eq. (fTTj) tells us that Ipoi,mm is given by P^ 



em. mm 



From Eqs. (fT4|) and p5|) . we have 



^pol.mm — ^em^mm ^mm 

^^ {-^ahs^opt/ ^ ) ^n 
C^ ^ahSjOpt ■ 



em, opt ^mm 



tmies 



(18) 



Here we use the fact r oc I^m ■ Note that the constant of 
proportionality does not affects the shape of the angular 
spectrum if grain properties do not vary much in halo. 
Therefore, as to the power spectrum q of Ipoi,mm, we 
can use that of Pabs^pUcaV- 



L-i Ot lpol,mm OC Cl of Pabs,opt- 



(19) 



Once we know the angular spectrum of Pabs .optical , we 
can estimate the angular spectrum of Ipoi,mm- 

Angular spectrum of starlight polarization, 
^abs, opt (l,b), is already available. Fosalba et al. (2002) 
obtained Ci ~ l~-^-^ for starlight polarization. The 
stars used for the calculation are at different distances 
from the observer and most of the stars are nearby 

Grain alignment theory based on radiative torques (Dolginov 
Sz Mytrophanov 1976, Draine & Weingartner 1996, 1997; Weingart- 
ner & Draine 2003; Lazarian & Hoang 2007ab; Hoang & Lazarian 
2008, 2009ab; see also Lazarian 2007, 2009 for reviews) predicts 
that grains starting with a particular size, which is ~ 5 X 10~® cm 
for the typical interstellar radiation field, get aligned. If the grain 
size distribution varies from one place to another, this fraction of 
aligned grains will also vary. However, for the diffuse ISM, the grain 
distribution does not vary much (Weingartner & Draine 2001). 
Therefore, the effect of variations of the degree of grain alignment 
due the variations of the grain size is of secondary importance. 
Much more important is the effect of incomplete alignment arising 
for some grains due to their shape and the direction of the radi- 
ation and magnetic field predicted in Lazarian & Hoang 2007a. 
The consequences of this effect for polarization maps require more 
studies (see Padoan et al. 2009). 

^^ We can make use of Right panel of Fig. [S] reversely. In the 
future, when we can accurately measure polarized emission from 
thermal dust in FIR or mm wavelengths, we can obtain the values 
of [{Cm,ax-Cmin)/iCmax +Cmin)]mm- This result Combined with 

the values of [{Cmax — Cmin)l(Cm,ax "t" Cmin)\o-ptical ™ Optical 

wavelengths can be used to find average angle between magnetic 
field and the plane of the sky. That is, when we know the ratio 
[...]mm/[...]optica!i we Can use Right panel of Fig. [5] to find the 
angle between magnetic field and the plane of the sky. 



stars. The sampled stars are mostly in the Galactic 
disk. CL02 reproduced the observed angular spectrum 
numerically using mixture of stars with a realistic 
distance distribution. CL02 also showed that the slope 
becomes shallower when only stars with a large fixed 
distance are used for the calculation. Therefore, it is 
clear that distance, or dust column density, to the stars 
is an important factor that determines the slope. We 
expect that, if we consider only the nearby stars with 
a fixed distance, the slope will be steeper. This means 
that, if wc consider stars in the Galactic halo, the slope 
will be steeper. 

The method described above requires measurements of 
polarization from many distant stars in the Galactic halo. 
Unfortunately, the number of stars outside the Galactic 
disk that can be used for this purpose are no more than a 
few thousands (Hciles 2000; see also discussions in Page 
et al. 2007; Dunkley et al. 2008). When more obser- 
vations are available, accurate estimation of Ipoi,mm{l, b) 
(and Ci of Pabs,opticai) will be possible. 

4.4. Model calculations for starlight polarization 

We expect the fluctuations of the starlight polariza- 
tion to arise primary from the fluctuations of magnetic 
fieldtH- The latter are expected to have the spectral 
index close to the Kolmogorov one (see §3). 

As we discussed in the previous subsection, the angular 
spectrum, Ci, of Pabs,opticai for the Galactic halo will be 
different from the observed l~^-^ spectrum for mixture of 
stars with different distances in the Galactic disk. How- 
ever, it is not clear exactly how the former is different 
from the latter. To deal with this problem we use nu- 
merical simulations again. We first generate two sets of 
magnetic field on a two-dimensional plane (8192 x 8192 
grid points), using Kolmogorov three-dimensional spec- 
trePn. Since we need Pabs. optical for stars well above the 
Galactic disk, we assume the distance to stars is fixed in 
each model. We consider 3 models: 

1. Case 1, Nearby stars in a homogeneous turbu- 
lent medium: We generate three (i.e. x,y, and z) 
components of magnetic field on a two-dimensional 
plane (8192 x 8192 grid points representing 400pc x 
400pc), using the following Kolmogorov three- 
dimensional spectrum: £'3_D(fc) ex k^^^'^ if fc > fco: 
where fco ~ 1/100 pc. (The outer scale of turbu- 
lence is lOOpc.) We assume the volume density of 
dust is homogeneous. All stars are at a fixed dis- 
tance of lOOpc from the observer. 

2. Case 2, Distant stars in a homogeneous turbulent 
medium: We generate three (i.e. x,y, and z) compo- 
nents of magnetic field on a two-dimensional plane 
(8192 X 8192 grid points representing 4 kpc x 4 kpc), 

^^ In addition, fiuctuations arising from the variations of the 
degree of grain alignment (see Lazarian 2007) are expected. 

^■^ Consider a 3-dimensional magnetic field with a 3D spec- 
trum E^£){kx,ky,kz) (ex: fc"™ with m = 11/3 for Kolmogorov 
turbulence). The spectrum of the magnetic field on a two- 
dimensional sub-plane (e.g. z = plane) is E^—q piane{kx,ky) oc 
J^ dkz E3i:>{kx,ky,kz), which we use to generate two sets of 
magnetic field on a two-dimensional plane. Note that, although 
this spectrum does not follow a power law near the outer scale of 
turbulence, it is close to k~"^~^^ for large values of k. 
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using the following Kolmogorov three-dimensional 
spectrum: E^jjik) oc fc^^^/^ if fc > fcoj where 
/co ~ 1/100 pc. Other setups are the same as those 
of Case 1, but the distance to the stars is 2kpc. 

3. Case 3, Stars in a stratified medium: We use 
the magnetic field generated in Case 1. The vol- 
ume density of dust shows a sech?'{z) decrease: 
p[r) = 4po/[exp(r/ro) -I- exp(— r/ro)]^. We assume 
spherical geometry and tq = lOOpc. The stars are 
at r = 200pc from the observer. The outer scale of 
turbulence is lOOpc. 

We assume that dust grains are oblate spheroids. In 
the presence of a magnetic field, some grains (especially 
large grains) are aligned with the magnetic field (see 
Lazarian 2007 for a review). Therefore, cross sections 
parallel to and perpendicular to the magnetic field are 
different. We assume that parallel cross section is ^30% 
smaller than the perpendicular one. We use the following 
equations to follow changes of Stokes parameters along 
the path: 

r^dI/ds^-5 + AaQ/I, (20) 

d{Q/I)/ds = A(j-Acr{Q/lf, (21) 

d{U/I)/ds = -Aa{Q/I){U/I) (22) 

(see Martin 1974 for original equations; see also Dolginov, 
Gnedin, & Silantev 1996), where 6 = {ai + 0-2), Aa = 
(cTi - 0-2), and 

2(Ti=f7i, (23) 

2(72 = cr± - (cr± - (T[|) cos 7 (24) 

(Lee & Draine 1985). Here a± and CT|| are the extinction 
coefficients and 7 is the angle between the magnetic field 
and the plane of the sky. After we get the final values of 
Stokes parameters, we calculate the degree of polariza- 
tion (■y/Q^ + C/^//) and, then, the second-order angular 
structure function of the degree of polarization. 

We show the result in Fig. [5] When all the stars are at 
the distance of lOOpc (Case 1), the spectrum is consistent 
with the Kolmogorov one for small 9. The result for 
the stratified medium (Case 3) also shows a spectrum 
compatible with the Kolmogorov one for small 9. When 
stars are far away (Case 2), the qualitative behavior is 
similar. However, if we measure average slope between 
9 — 0.2° and 20°, the result is different: the slope for 
Case 2 is substantially shallower. Note that 9 = 0.2° 
and 20° correspond to I = 1000 and 10, respectively. 
This means that, when we have distant stars only, the 
angular spectrum will be shallower than the Kolmogorov 
one. When we have a mixture of distant and nearby 
stars, we will have an angular spectrum that is steeper 
than the case of distant stars but shallower than the case 
of nearby stars, which implies the spectrum is shallower 
than Kolmogorov one. Therefore, it is not surprising 
that Fosalba et al. (2002) obtained a shallow spectrum 
of ~ C; oc Z~^-^ for a mixture of nearby and distant stars 
mostly in the Galactic disk. Flattening of spectrum (i.e. 
Ci ex ;-" with a w 1.3 - 1.4) for polarized FIR dust 
thermal emission is also observed in Prunet et al. (1998; 
see also Prunet & Lazarian 1999). 

Note that the spectrum of the emission polarization for 
very large values of multipole I is expected to be steeper 



than the spectrum of starlight polarization measured at 
small values of / For instance, from our model calcula- 
tions predict that we should see the Kolmogorov spec- 
trum of polarized emission, rather than ~ C; oc l~^'^. 
This difference must be taken into account if starlight po- 
larization is used to filter the polarized microwave emis- 
sion arising from dust. In fact, to do filtering one should 
either make a model of the magnetic field distribution 
based on the extended samples of stars throughout the 
Galactic volume or use only distant stars to get the spec- 
trum of polarization similar to that expected at the mi- 
crowave range. 

There could be systematic errors in transforming from 
the starlight polarization spectrum to that of emission. 
One of such errors may arise from the variation of the 
direction of the mean magnetic field direction and the 
line of sight. We expect that, while the starlight po- 
larization spectrum is relatively insensitive to the mean 
magnetic field direction, the emission spectrum shows a 
stronger dependence on it. Our preliminary calculations 
show that such a systematic error is small. We will pur- 
sue this possibility in the future. 

4.5. Comparison with the CMB polarization 

On the right panel of Fig. |51 we plot angular power 
spectrum of star light polarization. As we mentioned 
earlier, the angular spectrum of the degree of starlight 
polarization should be similar to that of polarized ther- 
mal dust emission (see Eqs. [18] and [T9]). 

To obtain angular spectra, we use a Gauss-Legendre 
quadrature integration method as described in Szapudi 
et al. (2001). To be specific, we first generate mag- 
netic fields from the 3 models we considered in the pre- 
vious subsection. Then, we calculate angular correla- 
tion functions, K{cos9). Finally, we obtain the angular 
spectra using Eq. ([T0|) . Since Ci obtained in this way is 
very noisy, we plot Ci averaged over the multipole range 
(//1-09, 1.090- We do not show Ci for I > 1000 because it 
is too noisy even with the averaging process. The second- 
order structure function on the left panel of Fig.[6]implies 
that l{l + l)Ci oc r5/3 foj. I > 2000. The straight dashed 
line for I > 1000 reflects this implication. 

We normalize the spectra using the condition ^2 (' + 
l)C;/27r — 3{fj,K^) for the acse of stars in the stratified 
medium. This normalization is based on the values given 
in Page et al. (2007; their Eq. [25] jQ. We assume that 
the observed band is W-band {v — 9AGHz). 

The plot shows that the slopes for nearby stars (thick 
solid line) and stars in the stratified medium are shal- 
lower than that of the Kolmogorov spectrum for / < 
1000. This result is consistent with that obtained with 
the angular structure function. Note that the case of 
distant stars has much flatter spectrum for I < 1000. 
We believe that our toy model for the stratified medium 
(Case 3) better represents the actual situation for po- 
larized emission from thermal dust in the Galactic halo. 
Therefore, we expect that the polarized thermal emission 

^■^ The purpose of this normalization is to match roughly our 
spectrum and that of Page et al. (2007) on large scales. Eq. (25) of 
Page et al. (2007) reads (l + l)Cf"=y2TT = 1.0{iy/65GHzfr^'^ - 
31-^-^ (fiK^) for 1/ = 94 GHz. When summed from I = 2 to I = 10, 
this gives ~ 2M{iiK^). We adopt J2l°{l + 1)Ci/2tt = 1,{iJ.K'^) for 
the case of stars in the stratified medium. 
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Fig. 6. — Model calculations for starlight polarization. Figures show the second-order angular structure functions (left) and the angular 
power spectra (right) for degree of polarization (right). Left: The second-order structure functions from model calculations. The structure 
function for nearby stars (thick solid line; d = lOOpc; Case 1 in the text) has a slope close to Kolmogorov one for 9 < 0.2° and slightly 
shallower spectrum for 0.2° < 6 < 10°. The case for density stratification (dashed line; scale height = lOOpc; Case 3 in the text) shows 
also a similar slope. However, the case for far-away stars (dotted line; d = Ikpc; Case 2 in the text) shows a substantially sh allow er slope 
for 9 > 0.1°. Right: The angular power spectra from model calculations. We take the shape of the spectrum from Section i|4.4l and the 
amplitude from Page at al. (2007) result. We assume that the observed frequency is 9AGHz. Nearby stars (thick solid line) or stars in 
stratified medium (dashed line) show a slope flatter than the Kolmogorov when / < 1000. The second-order structure function on the 
left panel implies that 1(1 + l)Ci oc l'^/^ for I > 1000. The straight dashed line for I > 1000 is obtained this way. Note that 9 = 0.1° 
corresponds to / ~ 180° /9° ~ 1800. We also show the polarized CMB 'EE' spectrum. 



from thermal dust in high-latitude Galactic halo has a 
spectrum slightly shallower than the Kolmogorov spec- 
trum for I < 1000. Our calculation do not tell us about 
the slopes for I > 1000. However, judging from the be- 
havior of the angular structure function for 9 < 0.1°, we 
expect that Ci oc Z"^^/'^ for / > 1000 (see the straight 
dashed line for I > 1000 on the right panel of Fig. [HI see 
also discussions in %A\i . 

We also show the polarized CMB 'EE' spectrum 
in the Figure (data from CMBFAST online tool at 
http://lambda.gsfc.nasa.gov/). The Figure shows that 
the EE spectrum dominates polarized thermal dust emis- 
sion from high-latitude Galactic halo for I > 100. The EE 
spectrum is expected to be sub-dominant when / > 5000. 

5. SPATIAL FILTERING OF FOREGROUNDS 

A unifying theme of this study and that in CL02 is 
that the principal source of the foreground fluctuations is 
related to the MHD turbulence in the Galactic interstel- 
lar medium. Our analysis of the available observational 
data in both publications supports this conclusion. In 
this section we discuss how this insight into the origin 
of the foreground fluctuations can be used to remove the 
foregrounds. 

5.1. Statistical properties of foregrounds 

Removal of Galactic foregrounds has always been a ma- 
jor concern for CMB studies. The challenge is only going 
to increase substantially now, when CMB polarization 
studies are attempted. 

The knowledge that foregrounds are not an arbitrary 
noise, but have well defined statistical properties in terms 
of their spatial power spectra is an important additional 
information that can be utilized to evaluate and eventu- 
ally eliminate the foreground contribution. 

Utilizing the information about underlying turbulence 
power spectrum is not straightforward, however. Our 
study shows that the observed power spectrum may de- 



pend on geometry of the emitting volume. Therefore, the 
detailed modeling of the foreground fluctuations should 
involve accounting for the geometry of the emitting vol- 
ume. 

The latter point stresses the synergy of Galactic fore- 
ground and CMB studies. Indeed, our fitting of the 
power spectra in Fig.[3]shows that on the basis of its vari- 
ations we may distinguish between different models of the 
emitting turbulent volume. As soon as this achieved, one 
can predict, for instance the level of fiuctuations that are 
expected from the foreground at the scales smaller than 
those studied. 

While the previous statements are true in general 
terms, a number of special cases in which simpler analysis 
is applicable are available. For instance, a simplification 
that is expected at higher resolutions that are currently 
available, is that at sufficiently small scales the statistics 
should get independent of the large-scale distributions of 
the emitting matter. Moreover, simple power laws are 
expected and observed (see Figure [T]) for the foregrounds 
at high Galactic latitudes. Therefore modeling of the dis- 
tribution of the Galactic emission and the filtering of it 
using the approach in §3 are not necessarily interlocked 
problems. 

5.2. Examples of filtering 

We remind the reader that the approach in §3 requires 
first separating a particular component of a foreground 
over a range of scaleJ^. e.g. using the traditional tech- 
nique of frequency templates, and then extending the 
spatial scaling of the foreground's Cf to higher I. Con- 
sider a few examples of utilizing this approach. 

For instance, high resolution measurements of the 

^^ If we know the foreground at a single spatial scale, one can 
use the a priori knowledge of the expected spatial scaling of the 
foreground spectrum. For instance, if underlying fiuctuations are 
related to magnetic field they are likely to have the spectrum close 
to the Kolmogorov one. 
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South Pole Telescope (SPT) (see Lueker et al. 2009) and 
the Atacama Cosmology Telescope (ACT) (see Fowler 
et al. 2010) will provide measurements at high I, but 
will have limited frequency coverage to remove the fore- 
grounds. Therefore, the spatial extrapolation of Cf 
obtained with other low resolution experiments, which, 
however, provide good frequency coverage required for 
the cf identification may be advantageous. 

Prior to the release of the Planck data, the exten- 
sion of cf obtained, for instance, with the WMAP to 
higher I may be useful for the foreground filtering for 
the suborbital missions. After the release of the Planck 
data, the studies of the advocated approach can still be 
useful. Consider, for instance. Fig. [6] If Planck mea- 
sures the spatial spectrum up to Z = 2000, then for a 
higher resolution balloon mission one can evaluate the 
level of foreground contamination by extrapolating the 
expected foreground spectrum. Moreover, the study of 
the B-modcs would require new experiments with higher 
sensitivity and the extrapolation procedure can be again 
useful. 

5.3. Demonstration of filtering technique 1 

In this subsection, we demonstrate how the filtering 
process works. For simplicity we use a Cartesian coor- 
dinate system. For demonstration purposes, we consider 
the CMB EE spectrum (see the thin-solid line in Fig.[6]3) 
and Galactic dust foreground for the case of stratified 
medium (see the dashed line in Fig. [6}d). 

Suppose that low-resolution all-sky maps are already 
available for many frequency channels. Let us assume the 
resolution of the maps is A6' ~ 21', which is similar to 
the WMAP resolution. Since there are many channels, 
one may use the usual filtering techniques to separate 
the CMB and the foreground signals. However, since we 
already know both the CMB and the foreground spectra 
in advance in this example, we do not follow the usual 
filtering techniques. Instead, we just assume that we 
already know the CMB and the foreground spectra for 
I < 500. 

The 'all-sky' map in the upper- left panel of Figure [7] 
represents this low-resolution map. The map is defined 
on a grid of 1024 x 1024 points. The angular resolution 
of the map is A0 ~ 21'. We generate the map from 
a much higher-resolution all-sky map (see upper-middle 
panel), which has a dimension of 8196 x 8196 in this 
example. When we obtain the low-resolution all-sky map 
from the much higher resolution all-sky map, we apply 
a circular top-hat beam pattern with radius ~ 10.5' to 
mimic actual observation. We calculate the spectra of the 
CMB and the foreground by direct Fourier transform of 
the low-resolution map data. 

Then, suppose that we perform a high-resolution bal- 
loon experiment that can cover only part of the sky. In 
such balloon experiments, frequency channels are usually 
limited and removing foregrounds is a challenging task. 

Here we show that, if we know the foreground spec- 
trum, we can easily obtain the CMB spectrum from the 
balloon data. The upper-right panel of Figure [7] repre- 
sents the high-resolution balloon data. The size of the 
data is 128 x 128. The angular resolution is /S.9 ^ 5.2' 
and the map covers 11° x 11° in the sky. We obtain the 
map by skipping every other point in each direction on 
low-right corner of the original map with 8192 x 8192 



points. 

Our goal is to obtain the angular spectrum of the CMB 
signal using the balloon data. We first make the balloon 
data periodic by proper reflections and translations and 
obtain a periodic data on a grid of 4096 x 4096 points. 
Then we multiply the data by a Gaussian profile of width 
~ 11°. The center of the Gaussian profile should locate 
near the center of the newly constructed 4096 x 4096 pixel 
data and coincide with the center of a 128 x 128 pixel 
original data. Then, we perform Fourier transform of the 
resulting data of size 4096 x 4096. This way, we obtain 
an angular spectrum of the total (i.e. CMB-|-foreground) 
fluctuations. 

Now, it is time to derive an angular spectrum of the 
foreground. We first take the angular spectrum from 
the low-resolution map, which is already available (see 
lower- left panel of Figure [7]) . The dashed curve is the 
foreground spectrum. The foreground spectrum is not 
defined for I > 500. In principle, we know the foreground 
spectrum when we know geometry and turbulence spec- 
trum. Here, we simply assume that the spectrum for 
/ > 500 is Kolmogorov. This way, we can construct the 
foreground spectrum for all values of I. 

The remaining task is just subtraction: when we sub- 
tract the foreground spectrum from that of the total fiuc- 
tuations, we can get the CMB spectrum. The solid curve 
in Figure [8] is the CMB spectrum obtained this way. The 
solid curve show a very good agreements with the original 
CMB spectrum used for generating the original all-sky 
data on a grid of 8192 x 8192 points (see middle panels). 
The r.m.s. relative percentage error (100AC;/Cp*^^) 
for 100 < I < 2000 is 24%, where AC/ is the differ- 
ence between the estimated CMB spectrum {(j^sumated .^ 

solid curve) and the true CMB spectrum (Cf^^^; dotted 
curve) . 

Note that the foreground spectrum depends on geom- 
etry of the emitting regions and underlying turbulence 
spectrum. In earlier sections, we discussed geometry of 
light-emitting regions. The turbulence spectrum may 
have some uncertainties (see Introduction) . Kolmogorov 
spectrum may be good for a number of cases. However, it 
is possible that the spectrum deviated from Kolmogorov 
one in some regions, such as molecular clouds. However, 
such deviation will not affect our filtering process much 
if the resolution of balloon data is not far beyond that of 
the low-resolution all-sky map. 

5.4. Demonstration of the filtering technique 2 

In the previous subsection, we demonstrated a recon- 
struction of CMB spectrum. The result is very promis- 
ing because the dust foreground model we used is not far 
from a realistic one. Note that we normalized the am- 
plitude of the dust foreground spectrum using the value 
quoted in an earlier work (Page et al. 2007). 

However, there are a couple of issues regarding to the 
previous demonstration. First, in the previous demon- 
stration, we assumed that we know the exact power-law 
index of the foreground spectrum. In general, we may 
not know the exact power-law index of the foregrounds. 
When this is the case, we need to extrapolate the fore- 
ground spectrum to larger values of multipoles and figure 
out to what extent we can safely extrapolate. Second, it 
is necessary to test how well our method works when the 
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Fig. 7. — Explanation of our filtering process. We use a Cartesian coordinate system for simplicity and assume there is only one 
foreground. Left panels: Low angular resolution all-sky map (1024 X 1024 pixels; A9 ~ 21'). We generate the low resolution map by 
applying a circular top-hat filter on the data shown in the middle panels. The radius of the filter is 4 grid units, which is equivalent to 
~ 10.5' The size of the resulting data file is 1024 X 1024, the angular resolution of which corresponds to A6 ~ 21'. We obtain the spectra in 
the lower-middle panel by direct Fourier transform. Middle panels: All-sky map of the CMB and the foreground signals and their angular 
spectra. We generate the map in the upper-middle panel using the spectra shown in lower-middle panel. The actual resolution of the map 
we use for the calculation is 8192 X 8192. Right panels: High angular resolution partial-sky map (128 X 128 pixels; A6 ~ 5.2'). This data file 
is from the low-right corner of the original 8192 X 8192 data. We obtain this map by skipping every other point in each direction, so that 
the angular resolution of this data {A8 ~ 5.2') is twice worse than the original 8192 X 8192 data, but 4 times better than the low-resolution 
data in middle panels. On the original data file (i.e. on the 8192 X 8192 data), the region spans 256 X 256 grid points. Therefore, the size 
of the partial-sky map is 128 X 128. The angular size of map in the upper-right panel is 11° X 11°. The maps shown in this figure are 
degraded to reduce file size. 
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Fig. 8. — Angular spectrum of CMB obtained from the filtering 
process. Top: The solid line is the spectrum obtained from the 
filtering process and the dotted line is the original CMB spectrum. 
Note that, since the high-resolution data span only 11° X 11°, the 
estimated angular spectrum for I < 100 is meaningless. Bottom: 
The solid line represents relative percentage error between the es- 
timated CMB (solid curve) and the true CMB (dotted curve) spec- 
tra. 

CMB and the foreground spectra are comparable. 



In this subsection, we develop a method of deriving 
a most probable extrapolation of the foreground spec- 
trum. We also test the performance of our technique for 
the case the CMB and the foreground spectra are com- 
parable. For this purpose, we scale up the amplitude of 
the foreground spectrum by a factor of 10000 

As in the previous subsection, we assume that multi- 
channel observations are available for I < 500. Therefore, 
we have foreground observation for I < 500. However, 
unlike the previous subsection, we assume that observa- 
tions are available only for selected values of multipole 
I. We mark 4 such I values on Fig. [HI We use the same 
foreground spectrum as in the previous subsection to gen- 
erate the 4 observed points. We also show the la obser- 
vation error bars for them, which are arbitrarily assigned 
in this example. 

Using the 4 observed points for foreground, we find the 
most probable extrapolation for I > 500. It is tempting 
to use the linear least-square fit in a log-log plane, in 

1^ The motivation we take the factor of 1000 is that the CMB 
BB spectrum by weak gravitational lensing is about 100-1000 times 
smaller than the CMB EE spectrum for I < 1000 (e.g. Zaldarriaga 
& Seljak 1998). 
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Fig. 9. — Second demonstration of the spatial filtering technique. For demonstration purposes, the foreground signal has been enhanced 
by 1000 times. Lejt: Most-likely foreground spectrum (solid line) and its la uncertainties (shaded region). We obtain the spectrum and 
the uncertainties using the observed data. In this example, we mark the 4 observed point with error bars. See the text for details. Right: 
CMB spectrum obtained from the spatial filtering technique. It shows that the spatial filtering technique works for I > 550, where the 
CMB and the foreground spectra have similar values. This result implies that the technique will work when the CMB spectrum is slightly 
larger than the foreground one. 



which X-axis is a logarithm of I and y-axis is that of 
l{l + l)Ci/2-K. But, it is not easy to estimate uncertainties 
in this case. 

Therefore, we adopt a shghtly different method. We 
work on the log-log plane. In this subsection, we use the 
following conventions: x = logl, y = \og[l{l + 1)Ci/2tt] 
and (xi. Hi) {i — 1, ..., N, where A^ = 4 in this example) 
denotes an observed point. For a given x, we want to 
find the most probable value and la uncertainty of y. 
In order to find them, we follow the following procedure. 
First, we select an arbitrary y and consider all possible 
lines that pass through {x,y). Let the slope of such a 
line be a and the y-intercept be b. Then, we calculate 
the following probability for the chosen y: 



P{y) = exp 



N 

E 



2a? 



(25) 



where z^ = axi + b and ct^ is the la observation error 
for i. Second, we repeat similar calculations for differ- 
ent values of y. Third, we find the value of y where 
Ptotiv) = HaiiUnesPiv) ^^^ ^he maximum. Fourth, we 
find la of the distribution Ptot{y)- Fifth, we repeat sim- 
ilar calculations for different values of x. 

We plot the result of this procedure in Fig. ^a). The 
solid line denotes the most probable values of 'y' and 
the shaded region represents the la uncertainty. The la 
uncertainty gets larger as I gets larger when I > 500. 
Therefore, it may not be a good idea to extrapolate for 
^s an order of magnitude larger than 500 which is the 
maximum I of the low-resolution multi-channel maps in 
our current example. However, it seems to be OK to 
extrapolate for Vs a few times larger than 500 in our 
current example. We can reduce the la uncertainty by 
reducing the observation error bars, which may enable 
us to extrapolate further for larger Vs. 

As in the previous subsection, we subtract the most- 
likely foreground spectrum from the observed (i.e. 
CMB-f-foreground) spectrum. We plot the result in 
Fig.Hb). 

Fig. m^b) shows that our technique recovers the CMB 
spectrum quite well for I > 550. Note that Fig.[9lja) tells 



us that the CMB and the foreground signals are almost 
same at I ^ 550. This example implies that our tech- 
nique works when the CMB spectrum is slightly larger 
than the foreground spectrum. It is not likely that our 
technique in the current form works if the foreground 
spectrum is larger than the CMB one. However, even in 
this case, it is possible that our technique can help to 
reduce observational uncertainties when combined with 
conventional multi-channel methods. 

The spatial filtering that we advocate here may be used 
as a part of a more general filtration procedure which 
uses both spatial and frequency information. Indeed, 
it is well recognized that the studies of tensor B-modes 
provide an excessively severe challenge to the precision 
of the removal of foregrounds. In this situation, it is 
important to reduce the errors of the determination of 
the foreground signal. The customarily used frequency 
templates provide filtering which is limited by both sys- 
tematic and measurement errors. In this situation, any 
additional information helping to decrease the errors is 
highly valuable. The spatial power spectrum of the fluc- 
tuations that we discussed in the paper does provide 
such an information. This approach can be illustrated 
with Figure 9, which currently has only points for low /. 
However, it is already can be seen that the initial error 
bars of the points can be decreased if we require that the 
points correspond to a power-law spectrum. If we had 
more points at higher I it is evident that the uncertainties 
could be further reduced. In other words, the constraint 
that the foreground fluctuations follow power law allows 
us to partially remove foregrounds at scales at which no 
foreground templates based on multi-frequency measure- 
ments are available and, at the same time, increase the 
precision of the foreground removal if the templates are 
available. 

6. DISCUSSION 

6.1. Our approach 

The major claim in this paper is that the statistical 
regularity of the foreground fluctuations enables one to 
extend their spectra from the scales where observations 
are available to the scales with no observations. In other 
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words, the spatial spectra of foregrounds are predictable. 
This, in its turn, makes possible spatial filtering of fore- 
grounds. 

The predictability of Galactic foreground fluctuations 
stems from the fact that they are due to ubiquitous 
Galactic turbulence. MHD turbulence is known to have 
well-defined statistical properties both in compressible 
and incompressible limits (see Cho & Lazarian 2005 for 
a review). Thus one expects to see power- law behavior, 
which in the case of magnetic field is expected to cor- 
respond to 3D power spectrum close to the Kolmogorov 
one. In the case of density, the spectrum is Kolmogorov 
for the subsonic turbulence, but gets shallower than the 
Kolmogorov one for supersonic turbulence (see ij4.ip . Ir- 
respectively of the underlying 3D spectrum being Kol- 
mogorov or not we claim that one can predict the entire 
spectrum of spatial 2D foreground fiuctuations, when the 
measurements over a limited range of scales are available. 
For this purpose we need to know the geometric proper- 
ties of the volume. Our expectations of the underlying 
spectrum to be in some instances, e.g. for magnetic field, 
to be close to Kolmogorov only help to increase the ac- 
curacy of our prediction of the 2D foreground spectrum. 

With the 2D predicted spectrum one then can use 
the procedure of filtering the foreground using Eq. ^ 
(see ij5.3l for demonstration of the procedure) . Alterna- 
tively, the good correspondence of the observed spatial 
spectrum of the foregrounds may serve as an additional 
proof of the accuracy of the foreground removal proce- 
dure. Needless to say, the information of the Galactic 
turbulence and the geometry of the emitting region that 
can be a by-product of the CMB research is of astrophys- 
ical significance. 

A note of warning is due, however. The procedure of 
statistical filtering that we demonstrated in the paper, 
as any other procedure of the foreground removal, is not 
ideal. The power-law approximation of Cf is definitely 
not exact. In this paper we have analyzed the causes for 
such deviations and provided the explanation for most 
notable features characterizing the change. More de- 
tailed modeling of Galactic turbulence is required. 

6.2. Our data analysis 

In this paper, we have discussed angular spectra of 
Galactic foregrounds. We have focused on synchrotron 
total intensity and polarized thermal dust emission. Our 
current study, as well as earlier studies (Chepurnov 1999; 
CL02), predicts that C; will reveal true 3D turbulence 
spectrum on small angular scales. 

Our model calculations that take into account stratifi- 
cation effects imply that 

1. 6 < a. few times 0.1° (or / > a few times 100) for 
synchrotron emission (see Fig. [3]), and 

2. 6 < 0.1° (or / > 1000) for polarized emission from 
thermal dust (see Fig. [6|). 

On larger angular scales, spectra are expected to be shal- 
lower. 

6.3. Results 

In this paper, we have analyzed Haslam 408MIIz map 
and a model dust emission map and compared the results 
with model calculations. We have found that 



1. The Haslam map for high Galactic latitude (6 > 
30°) can be explained by MHD turbulence in the 
Galactic halo. The measured second-order an- 
gular structure function is proportional to 9^'^, 
which corresponds to an angular spectrum of l~^'^. 
The high-order statistics for high Galactic latitude 
{b > 30°) is consistent with that of incompressible 
magnetohydrodynamic turbulence. Our model cal- 
culations show that a two-component model (see 
§3. 31 and Fig.[3l) can naturally explain the observed 
angular spectrum. The one-component model can 
also explain the observed slope. But, the slope 
of the spectrum shows a more abrupt change near 
Z-30. 

2. The model dust emission map may not have any- 
thing to do with turbulence on large angular scales. 
That is, we do not find signatures of turbulence in 
the map. 

3. Both maps show flat high-order structure functions 
for the Galactic plane. This kind of behavior is ex- 
pected when discrete structures dominate the map. 

We have described how we can obtain angular spec- 
trum of polarized emission from thermal dust in high 
Galactic latitude regions. Our model calculations show 
that starlight polarization arising from dust in high 
Galactic latitude regions will have a Kolmogorov spec- 
trum, Ci oc Z^^^/^, for I > 1000 and a shallower spectrum 
for / < 1000 (Fig. E}. We expect that polarized emission 
from the same dust also has a similar angular spectrum. 
That is, we expect that the angular spectrum of polar- 
ized emission from thermal dust is close to a Kolmogorov 
one for I > 1000. 

We have described a new technique of filtering CMB 
foregrounds. When we have 

1. low angular resolution full-sky measurements (such 
as WMAP data) and 

2. high angular resolution partial-sky measurements 
with limited frequency channels, 

we can use the technique to derive the CMB angular 
spectrum for the high- resolution data. In §5.31 and HbAl 
we have demonstrated the technique. 

6.4. Comparison with approaches in the literature 

The existing confusion in the literature include naive 
identification of the 2D spectra of foregrounds with the 
spectra of the underlying fluctuations. Therefore, for 
instance, from the fact that the spectral slope of Cf 
differ from the Kolmogorov one, the conclusion about 
the nature of the fiuctuations is made. 

In the paper above we have shown that the 2D spectra 
may have a spectral slope different from the underlying 
spectral slope of the turbulence. We showed that it is 
essential to take into account the non-trivial geometry 
of observations with the observer sampling turbulence 
along the diverging lines of sight and within the volume 
where the density of emitters changes. For the case of 
the synchrotron fluctuations we showed that the observed 
non-Kolmogorov value of the spectral index of Cf can 
be reconciled with the Kolmogorov- type turbulence. 
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A notable difference of our study compared with CL02 
is that in the latter study we tried to necessarily asso- 
ciate the spectra of foregrounds with the underlying Kol- 
mogorov or Goldreich-Sridhar (1995) turbulence with the 
spectral slope —11/3. The limitations of these approach 
get evident in view of establishing of the shallow spectral 
index of density fluctuations in supersonic MHD turbu- 
lence (Beresnyak et al. 2005). These fluctuations, ac- 
cording to our present study, can explain the observed 
spectra of Galactic dust emission. 

6.5. Limitations and extensions of our filtering 
procedure 

The procedure of spatial statistical filtering discussed 
in §2 is very simplified. It should provide satisfactory re- 
sults for a simple power law behavior of Cf (see Figure 
2). In more complex cases detailed modeling of the emit- 
ting volume and/or use of the filtering as a part of a more 
sophisticated foreground removal procedure is required. 
We partially addressed this issue in SJH Especially in 
N5.41 we described a method to estimate uncertainties 
stemming from a power-law modeling of the foreground 
spectrum. Nevertheless, more improvement is needed for 
the new technique. At the same time, we believe that our 
approach can be applied to the removal of foregrounds 
not only from the CMB data, but, for instance, from the 
high-z hydrogen statistics studies. 

7. SUMMARY 

In the paper above we have obtained the following re- 
sults: 

1. We provided additional evidence that the syn- 
chrotron Galactic emissivity is consistent with the halo -f 
disk model. Within this model we show that the spatial 
spectrum of the underlying 3D fiuctuations is consistent 
with the Kolmogorov one. 

2. Within our model we related the angular scale for 
the change of the power spectrum of the synchrotron 
fluctuations with the ratio of the injection scale to the 
thickness of the observed region in the direction of ob- 
servation. 

3. We explained the spectrum of dust foreground emis- 
sion as arising from the shallow spectrum of density 
fluctuations which characterizes supersonic MHD turbu- 
lence. 



4. We used numerical modeling of grain optical prop- 
erties to relate the polarization of starlight with the ex- 
pected sub-mm foreground polarization and outlined the 
ways of quantitative use of starlight polarization maps 
to study sub-mm polarized dust foreground. We evalu- 
ated the uncertainty of the evaluated sub-mm polariza- 
tion spectrum arising from the variations of the magnetic 
fleld direction in respect to the line of sight. 

5. We showed that for randomly chosen sample of stars 
the spectrum of the starlight polarization for the under- 
lying Kolmogorov turbulence depended on how the se- 
lected stars are distributed along the line of sight. For 
stratified model of Galactic dust we predicted the spec- 
trum of spatial fluctuations with the index approaching 
the Kolmogorov value of —11/3 for I > 1000. 

6. On the basis of our improved understanding of the 
self-similar nature of the underlying MHD turbulence we 
proposed and tested a procedure of spatial statistical re- 
moval of Galactic foregrounds based on extending of the 
Cf spectrum to higher I. 

7. We studied the higher order correlations of the 
Galactic dust and synchrotron foregrounds and reported 
substantial difference with the higher-order scalings (see 
Appendix D). The synchrotron scaling shows intcrmit- 
tencies similar to magnetic fields in incompressible fiuids, 
while dust scaling is closer to the intermittency demon- 
strated by highly compressible MHD turbulence. 
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APPENDIX 



SPECTRUM AND STRUCTURE FUNCTION OF DIFFUSE SYNCHROTRON EMISSION: A MODEL FOR HOMOGENEOUS 

TURBULENCE REVISITED 

Suppose that 3-dimensional (3D) MHD turbulence has a 3D spatial power spectrum of the form E^d oc fc~™, where 
k is the wavenumber. Note that in Kolmogorov turbulence m = 11/3. Then what will be the 2-dimensional angular 
spectrum, Ci, of the observed synchrotron total intensity? We cannot not directly observe the 3D spatial power 
spectrum. However, we can infer the 3D spectrum from the observed 2-dimensional (2D) angular spectrum, Ci, of the 
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synchrotron total intensity. Then let us find out how the 3D spectrum, E^oik), and the 2D angular spectrum, C/, are 
related. We mostly follow discussions in Cho & Lazarian (CL02). We also make use of an analytical insight obtained 
in Lazarian (1992, 1995ab) and numerical results obtained in CL02. Although we focus on synchrotron emission here, 
the discussion in this section can be applicable to any kind of emission from an optically thin medium. 

MHD turbulence and synchrotron emission 

For synchrotron radiation, cmissivity at a point r is given by e(r) oc n(e)|_Bj^|''', where n(e) is the electron number 
density, B± is the component of magnetic field perpendicular to the line of sight. The index 7 is approximately 2 for 
radio synchrotron frequencies (see Smoot 1999). If electrons are uniformly distributed over the scales of magnetic field 
inhomogencities, the spectrum of synchrotron intensity reflects the statistics of magnetic field. For small amplitude 
perturbations {5b/ B <i; 1; this is true for scales several times smaller than the outer scale of turbulence if we interpret 
B as local mean magnetic field strength and Sb as random fluctuating fleld in the local region), if Sb has a power- 
law behavior, the synchrotron emissivity will have the same power-law behavior (see Getmantsev 1959; Lazarian & 
Shutenkov 1990; Chepurnov 1999). Therefore, we expect that the angular spectrum of synchrotron intensity also 
reflects the spectrum of 3-dimensional MHD turbulence. 

When an observer is located inside a turbulent medium, the angular correlation function, hence the power spectrum, 
shows two asymptotic behaviors. When the angle is larger than a critical angle, we can show that the angular 
correlation shows a universal 9~^ scaling. On the contrary, when the angular separation is smaller than the critical 
angle, the angular correlation reflects statistics of turbulence. In this small angle limit, we can show that the angular 
power spectrum is very similar to that of turbulence. The critical angle is determined by the geometry. Let the outer 
scale of turbulence be L and the distance to the farthest eddies be dmax- Then the critical angle is 

e - L/d^ax- (Al) 

Small- angle limit in homogeneous turbulence 

When the angle between the lines of sight is small (i.e. 6 < L/ dmax), the angular spectrum Ci has the same slope as 
the 3-dimensional energy spectrum of turbulence. Lazarian & Shutenkov (1990) showed that if the 3D spatial spectrum 
of a variable follows a power law, E^oik) oc /c"™, then the 2-dimensional spectrum of the variable projected on the 
sky also follows the same power law, 

Ci ex r™ (A2) 

in the small 9 limit. For Kolmogorov turbulence {E^u oc k^^^'^), we expect 

Gocr"/3, iie<L/dmax- (A3) 

Note that I ~ T:/e. 

In some cases, when we have data with incomplete sky coverage, we need to infer C'l from the observation of the 
angular correlation function 

w{e) =< /(ei)/(e2) >, (A4) 

where /(e) is the intensity of synchrotron emission, ei and 62 are unit vectors along the lines of sight, 9 is the angle 
between ei and 62 , and the angle brackets denote average taken over the observed region. As we discuss in Appendix 
A, when the underlying 3D turbulence spectrum is oc A:~™ (e.g. m = 11/3 for Kolmogorov turbulence), the angular 
correlation function w{9) is given by 

w{e) oc < /2 > -const 9""-^, if 61 < L/d^ax- (A5) 

It is sometimes inconvenient to use the angular correlation function in practice to study turbulence statistics because 
of the constant < P >. 
A better quantity in small-angle limit would be the second-order angular structure function: 

D2{9) = <\I(ei)-I(e2)\^> (A6) 

= 2<P>-2w{9). (A7) 
Thus, in homogeneous turbulence with 3D spatial spectrum of E{k) oc A:"™, we have 

D2i9) oc 9""-^. (A8) 

When we measure the slope of the angular structure function, we can infer the slope of the 3D spatial power spectrum 
of turbulence. 

Large-angle limit in homogeneous turbulence 

In this limit, the angular correlation function is more useful than the structure function. Following Lazarian & 
Shutenkov (1990), we can show that the angular correlation function for 9 > L/dmax follows 

w{9)^ f fdridr2lC{\ri-r2\), 
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Fig. 10. — Two limits in homogeneous turbulence. Upper plot: Small limit {9 < L/dmax)- The fluctuations along the entire length 
of the lines of sight are correlated. Lower plot: Large 8 limit (8 > L/dmax)- Only points close to the observer are correlated. Note the 
definition of r and ■(/>. From CL02. 

1 r^ K^/^^"'^'^, ^-^ const .^^. 

where lC{r) is the 3D spatial correlation function and we change variables: (ri, r2) ^ (r, V'), which is clear from Fig. llOl 
We accounted for the Jacobian of which is r/ svaO. We can understand \/9 behavior qualitatively as follows. When 
the angle is large, points along of the lines-of-sight near the observer are still correlated. These points extend from the 
observer over the distance oc l/sin(6'/2). 
If we assume L/dmax < ^ 1, we can get the angular power spectrum Ci using Fourier transform: 

Ci^ j fw{0)e-'''d0,d0y 



de ejo{ie)w{e) oc z-\ (aio) 

where 9 = (0^ + OyY^"^ , Jo is the Bessel function, and we use w{9) oc 6~^. 

Expectations for homogeneous turbulence 

In summary, for homogeneous Kolmogorov turbulence (i.e. E{k) oc fc^^^"), we expect from equations (|A3p and 
((Mo)) that 



which means that the power index a of Ci ifo — 1 < a < —11/3. We expect the following scaling for the second-order 
angular structure function: 

^^ -^ °^ I constant if 6* > L/dmax- 

The critical angle 9cr ~ L/dmax depends on the size of the large turbulent eddies and on the length of the line of sight. 
If we assume that turbulence is homogeneous along the lines of sight and has L ~ 100 pc corresponding to a typical 
size of the supernova remnant, and that dmax ~ 1 kpc for synchrotron halo (see Smoot 1999), we get 9cr ^ 6°. 

THE SECOND-ORDER ANGULAR STRUCTURE FUNCTION IN THE SMALL ANGLE LIMIT 

In this Appendix, we discuss how underlying 3D statistics and observed 2D statistics are related. This discussion 
is useful when we infer 3D statistics from observed 2D statistics, or vice versa. Strictly speaking, discussion in the 
section is applicable to homogeneous and isotropic turbulence onlyH 

The angular correlation w{9) is given by the integral 

w{9)^ JdhJdl2lC{\h-h\), (Bl) 

^^ Note that point sources would result in a ~ 0. j m-, i i i- r 4.- / ot-> ^ 4- i- c 

, c, ^ and 2D angular correlation function (or 2D angular structure func- 

^« When turbulence is inhomogeneous or anisotropic, we may ^j^^) discussed in this section is also applicable to inhomogeneous 

not directly apply the r esults m this section. However, our numeri- cases 

cal calculations in i|3.3l shows that the relation between 3D statistics 
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where /C(r) is the 3 dimensional spatial correlation. Suppose that li is along x-axis, li = {h, 0), and I2 — l2{cos9, sin 9). 
Then, the correlation in the limit of small 9 is given by 



dh f dh JC U{h~l2COs9Y + lls\ 



w{9)^ dh dhJCiJih- h cos 9f + ll sin" 9 (B2) 



dh dhK. [J {ll-2hl2 cos 9 + 11] (B3) 



dl 



dh IC [Vih-hy + hhO^) . (B4) 



Suppose that the spatial correlation follows a power law: /C(r) oc const — r^ for r < L, where L is the outer-scale of 
turbulence. For Kolmogorov turbulence, p — 2/3. Then the derivative of 'w{9) is given by 



w(9) 



I dh I dh [{h-hf+hhe'^Y^'' \2hhO) (B5) 

du I dw [u;^ + (u2-w2)5)2/4jP/2-i^^2_^2)g,/2, (B6) 

where u = h + h and w ~ h ^ h- If P < 1, the integration diverges as 9 goes to zeroH Therefore, when p < 1, it 
suffices to perform the integration in the vicinity of Zi = ^2 or w = 0. Then we have 

^(^).. Z"^., Z"^.,, L,.2 , .,2a2/.lP/2-l ,2 



(x~ du dw [u;" + u"6i74j""' u''9/2 (B7) 

du {u9/2)P-^u^9/2 oc -9P, (B8) 

where we use J_^ dw/{w^ + A^)" = yl^^^" /-W2 '^^ sec^"^" 9. Therefore, for small 9 we have 

w{9) oc Ci - C29P+\ (B9) 

where Ci and C2 are constants. Comparing this equation with 

w{9)^C3~CiD2i9), (BIO) 

we get 

D2{9) oc 9P+\ (Bll) 

Analytic expressions for the relation between the angular structure function (-D2) and the spatial ID spectrum (E(k); 
in case of Kolmogorov, E{k) oc fc~^/^) can be found in the literature. For example, Lazarian (1995a; see also Lazarian 
& Shutcnkov 1990) derived the following expression: 

/■^/■^ d 
E{k)ock dr] — {Q{f])f]) JiikRr]) + IC5, (B12) 

Jo dri 

where C can be regarded as the outer scale of turbulence, R is the size of the system, Q{ri) ^ D'2{rj)rj, rj — sm9, Ji{x) 
is the Bessel function of the first order, and JCq is a small correction term. 

SPATIAL SPECTRUM OF EMISSIVITY 

The synchrotron emissivity is proportional to ~ n{e)B'^ oc B^, where n(e) is the high-energy electron number density. 
Suppose that magnetic field is roughly a Gaussian random variable. This may not be exactly true, but should be a 
good approximation. When a Gaussian random variable _B(r)(2_| follows a Kolmogorov spectrum 

n I n^l \|2 f if fc < fco ,„,, 

Eb,3d = \Bik)\ oc I (fc/fc^)-ii/3 if fc > k,, (CI) 

we can show that the 3D spectrum of B^{r) follows Eq. IH] (see, for example, Chepurnov 1999). The correlation of 
B^ir) and 3D energy spectrum of B'^{r) are related by 

Kb^Iv) =< i32(x)B2(x + r) >,^ J £1,2^^)6'''"' d% (C2) 

^^ When p = 1, the spatial correlation becomes IC(r) oc C — r, ,, ,■ ■ , , r i u 1 1,11 ti, 7-4 

, ^ ■ , J mi ,■ 01-^ J ■ n/7N three-dimensional spectrum oi turbulence is shallower than k *, 

where C is a constant. The corresponding 3D spectrum is E{k) oc ^^ ^^^^ ^^^^ ^ ^^ _ ^„_3_ ^^^^^ ^^ ^ ^m-3 ;, ^ constant. 

fc_ When the slope of the turbulence spectrum is steeper than Therefore, the condition of p < 1 is generally satisfied in turbulent 

k , the correlation function has the form /C(r) oc /Co — r re- medium 

gardless of the turbulence slope. On the other hand, when the 20 For simplicity, we assume B is a scalar. 
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Eb-M^) = |B2(fc)|2 oc j KB2{v)e-'''-'<fv, (C3) 

where < ... >x denotes an average over x. A Gaussian random variable satisfies 

< B2(x)b2(x + r) >=< B2(x) >< B2(x + r) > +2 < B(x)B(x + r) >^ , (C4) 

where the first term on the right is a constant. Therefore we can ignore the term in what follows. Fourier transform 
of both sides results in 

LHS = SB2,3B(k), (C5) 

RHS = 2 / < B(x)B(x + r) >2 e-'^^'d^v (C6) 

= 2 / rf^r/CsW/CsWe-*'^"- (C7) 

= 2J Sv j d'pj d'qEB,3Dip)EBMny^''^^'''^"' (C8) 

= 2 / rf^p/" d^'q EB,3D{p)EBMn)S{p + q-k) (C9) 

- 2 / d^p £;b,3Z5(p) SB.3D(k - p), (CIO) 

where S{k) is the Dirac (5-function. Therefore we have 

SB^3i)(k) ~ ^52^3^(0) « 2 /" d^k\EB.M^}\^ ~ constant (CU) 

for fc <C fcg. 

HIGH-ORDER STATISTICS 

While most of the paper is directly related to making use of the knowledge of the underlying spectra and/or two 
points correlations in order to remove foregrounds, the part dealing with higher-order statistics is not directly related to 
the foreground removal. Nevertheless, the correspondence of the intermittencies of foregrounds to those of turbulence 
provide another support for our understanding of the turbulent origin of foreground fluctuations. 

High-order structure functions are used for the study of intermittency, which refers to the non-uniform distribution 
of structures. Since CMB signals are close to Gaussian, one may think that they do not have strong intermittency. 
However, there is a report that intermittency of CMB signals deviate from Gaussianity (Bershadskii & Screenivasan 
2003). If intermittency of foreground signals are different from that of CMB signals, we can potentially use high-order 
structure functions to separate CMB and foreground signals. In addition, we can use high-order structure functions 
to improve our knowledge of foregrounds. 

The structure functions of order p for an observable / is defined by 

Sp(r)^<\I{x)-I{x + r)\P>, (Dl) 

where the angled brackets denote average over position x. For an observable defined in the plane of the sky, the angular 
structure function of order p is 

Dp{e)=<\I{e,)-I{e2)\P>. (D2) 

Traditionally, researchers use high-order structure functions of velocity to probe dissipation structures of turbulence. 
In fully developed hydrodynamic turbulence, the (longitudinal) velocity structure functions Sp =< ([v(x-l-r) — v(x)] • 
f)P >=< Sv^{r) > are expected to scale as r^p; 

Sp{r) cc r'^" . (D3) 

One of the key issues in this field is the functional form of the scaling exponents Cp- There are several models for (p. 
Roughly speaking, the dimensionality of the energy dissipation structures plays an important role. 
Assuming 1-dimensional worm-like dissipation structures. She & Leveque (1994) proposed a scaling relation 

Cp^^=p/9 + 2[l- (2/3)^/3] (D4) 

for incompressible hydrodynamic turbulence. Note that (^p is the scaling exponent of structure of order p. On the 
other hand, assuming 2-dimensional sheet-like dissipation structures, Miiller & Biskamp (2000) proposed the relation 

Cf^=p/9 + l- (1/3)^/3 (D5) 
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for incompressible magneto-hydrodynamic turbulencq^j. 

Recently, high-order structure functions of molecular line intensities have been also employed (Padoan et al. 2003; 
Gustafsson et al. 2006). In optically thin case, the molecular line intensities are proportional to the column density. 
Kowal, Lazarian, & Beresnyak (2007) studied scaling of higher moments of density fluctuations in MHD turbulence. 
Their numerical results show that, first of all, the scalings of higher moments in 3D can be obtained by studying the 
distribution of column densities. Then, they showed that the behavior of the scaling exponents for column density 
depends on sonic Mach number of turbulence. 

Bershadskii & Screenivasan (2003) calculated intermittency of CMB signals. Their result shows that the WMAP 
data follow She-Leveque scaling. Therefore, it is interesting to see if foreground signals show different scalings. 

We plot the scaling exponents of the Haslam 408MHz map in Fig. [TTJ The Haslam map shows a reasonable 
agreement with the MuUer-Biskamp MHD model. Note that Padoan et al. (2003) also obtained a similar result using 
^^CO emission from Perseus and Taurus. 

Unlike the Haslam map, the 94GHz dust map does not show agreement with MuUer-Biskamp model. The scaling 
exponents do not show strong dependence on the order p similar, which makes them similar to the results of the 
potential part of the supersonic fiow studied in Kowal & Lazarian (2010) and densities in supersonic super-Alfvenic 
flows studied in Kowal et al. (2007). One may expect that in the presence of self-gravity this is what is expected for 
higher order statistics of density, but further studies are required. 

Left and middle panels of Fig. [12] show that the slope is around 1 for high-order structure functions. This kind 
of behavior is expected when discrete structures dominate the map. However, it is not clear what kinds of discrete 
structures dominate^- Since the Haslam map and the dust map sample different types of the ISM, it is not so 
surprising that they show different scaling behaviors. 

For the Galactic disk, high-order structure functions of both maps show nearly flat structure functions. The structure 
functions show a nearly flat behavior down to 9 ~ 1°, which happens to be the actual angular resolution of the 
maps. Further studies of MHD turbulence in the presence of self-gravity should clarify the origin of such behavior. 
Alternatively, higher moments of structure functions can be affected by unresolved point sources. 
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Fig. 11. — High order statistics of the Haslam 408MHz map. Left: The second, third, 5th, and 8th-order structure functions for b > 30°. 
Middle: The scahng exponents for b > 30° seem to follow the MuUer-Biskamp MHD scaling. Note that we equate the scaling exponent of 
the observed third-order structure function and that of the MuUer-Biskamp model (see more discussion in Cho et al. 2003). We measure 
the slope between = 1° and 10°. Right: The second, third, 5th, and 8th-order structure functions for the Galactic disk (—2° < 6 < 2o). 
They are all flat. 



^^ Boldyrev (2002) obtained the same scaling relation for highly 
supersonic turbulence. However, since it is unlikely that turbulence 
in the Galactic halo is highly supersonic, we refer the scaling rela- 
tion to the "MuUer-Biskamp" scaling. 

^^ However, we can show that thin filamentary structures or 
point sources with uniform intensity and sharp boundary are not 
the dominant structures. Consider a circular cloud with a ra- 
dius A and a uniform intensity / centered at the origin. (For 
simplicity, let us consider a two-dimensional cloud in the two- 
dimensional Cartesian coordinate system. Let us assume inten- 
sity is defined on the Cartesian grid.) Then the structure func- 
tion is given by Dn(r) oc J^iY.^ \H^i) - I{^j)l" /Npair, where 
two points Xi and Xj are separated by the distance r and Npai,- 
is the total number of such pairs. When r ^ A, we can show 



that Ei Ej \H^i) ~ J{^j)\" « /"(27rr)(7rA2) and Npa^r oc (27rr). 
Therefore, we have Dn{r) oc /"(irA^) for r 2> A. That is, the 
structure functions show no dependence on r. We expect that 
structure functions for a thin uniform filament also show no de- 
pendence on r when r is larger than the width of the filament, 
because a filament can be viewed as a chain of circular clouds (or a 
chain of square-like clouds). The high-order structure functions for 
the dust emission map show ~ r^ power-law scaling for 6 > 1° (left 
panel of Fig. Illll . Thus neither thin uniform filaments nor point 
sources are dominant structures. However, filaments or circular 
clouds with smoothly extended boundaries can be the dominant 
structures. 
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Fig. 12. — High order statistics of the 94GHz dust emission map. We use dust emission intensity times sin 6 for all calculations. Left: 
The second, third, 5th, and 8th-order structure functions for b > 30° . Middle: The scaling exponents for b > 30° does not show signatures 
of turbulence. We measure the slope between 9 = 1° and 10°. Right: The second, third, 5th, and 8th-ordcr structure functions for the 
Galactic disk (-2° < 6 < 2o). They are all flat. 



